text:
# The decay of homogeneous anisotropic turbulence

## I Introduction

Decaying turbulence has attracted the attention of various communities and is often considered in experimental, numerical and theoretical investigations . It is in fact quite common that even experiments aimed at studying stationary properties of turbulence involve processes of decay. Important examples are provided by a turbulent flow behind a grid (see  and references therein) or the turbulent flow created at the sudden stop of a grid periodically oscillating within a bounded box . In the former case, turbulence is slowly decaying going farther and farther away from the grid and its characteristic scale becomes larger and larger (see  for a thorough experimental investigation). Whenever there is sufficient separation between the grid-size $`L_{in}`$ and the scale of the tunnel or the tank $`L_0L_{in}`$, a series of interesting phenomenological predictions can be derived. For example, the decay of the two-point velocity correlation function, for both isotropic and anisotropic flows, can be obtained under the so-called self-preservation hypothesis (see  chapter XVI). That posits the existence of rescaling functions allowing to relate correlation functions at different spatial and temporal scales. By inserting this assumption into the equations of motion, asymptotic results can be obtained both for the final viscosity-dominated regime and for the intermediate asymptotics when nonlinear effects still play an important role.
The status of the self-preservation hypothesis and the properties of energy decay in unbounded flows are still controversial . Systematic results on related problems have been established recently, e.g. for non-linear models of Navier-Stokes equations as Burgers’ equation, see e.g. , and for stochastic models of linear passive advection , both in unbounded  and bounded domains .
Here, we investigate the decay of three-dimensional homogeneous and anisotropic turbulence by direct numerical simulations of the Navier-Stokes equations in a periodic box. Previous numerical studies have been limited to either homogeneous and isotropic turbulence  or to shell models of the energy cascade .
The initial conditions are taken from the stationary ensemble of a turbulent flow forced by a strongly anisotropic input . The correlation lengthscale of the initial velocity field $`L_{in}`$ is of the order of the size of the box $`L_0L_{in}`$.
In the first part of this paper, we shall try to answer the following questions about the intermediate asymptotic regime of nonlinear decay: How do global quantities, such as single-point velocity and vorticity correlations, decay ? What is the effect of the outer boundary on the decay laws ? Do those quantities keep track of the initial anisotropy ? As for the statistics of velocity differences within the inertial range of scales, is there a recovery of isotropy at large times ? If so, do strong fluctuations get isotropic at a faster/slower rate with respect to those of average intensity ? Do isotropic and anisotropic fluctuations decay self-similarly ? If not, do strong fluctuations decay slower or faster than typical ones ?
In the second part we study the early stages of the decay, with the aim of establishing a link between the small-scale velocity statistics in this phase and in the forced case. That will allow us to argue in favor of an “exponents only” universality scenario, for forced hydrodynamic turbulence.

## II Numerical setup

### A The initial conditions

The initial conditions are taken from the stationary ensemble of a forced random Kolmogorov flow . For sake of completeness, we recall here some of the statistical properties of this forced turbulent flow. We consider the solutions of the Navier-Stokes equations for an incompressible velocity field $`𝒗`$.

$$_t𝒗+(𝒗)𝒗=p+\nu \mathrm{\Delta }𝒗,$$
(1)
in a three-dimensional periodic domain. To maintain a statistically stationary state Eq. (1) had to be supplemented by an input term $`𝒇`$ acting at large scales. This force was strongly anisotropic: $`𝒇=(0,0,f_z(x))`$ with $`f_z(x)=F_1\mathrm{cos}[2\pi x/L_x+\varphi _1(t)]+F_2\mathrm{cos}[4\pi x/L_x+\varphi _2(t)]`$, constant amplitudes $`F_{1,2}`$ and independent, uniformly distributed, $`\delta `$-correlated in time random phases $`\varphi _{1,2}(t)`$. This choice ensured the statistical homogeneity of the forcing and thus of the velocity field. We simulated the forced random Kolmogorov flow at resolution $`256^3`$ for time spans up to $`70`$ eddy turnover times . The viscous term was replaced by a second-order hyperviscous term $`\nu \mathrm{\Delta }^2𝒗`$. We stored $`40`$ statistically independent configurations that here serve as initial conditions for the decaying runs.

### B Decaying runs

As turbulence decays, the effective Reynolds number $`Re=v_{rms}L_0/\nu `$ decreases, while the viscous characteristic scale and time increase. To speed up the numerical time-marching, it is then convenient to use an adaptative scheme. We calculate periodically the smallest eddy-turnover time from the energy spectrum and set the time step as $`1/100`$ thereof. The whole velocity-field configuration is then dumped for offline analysis at fixed multiples $`\{0,1,10,10^2,10^3,10^4,10^5,10^6\}\tau _0`$ of the initial large-scale eddy turnover time $`\tau _0=L_0/v_{rms}^{t=0}`$. In Fig. 1 we show a two dimensional section in the plane $`x`$-$`z`$ of the velocity components $`v_z`$ and $`v_x`$.

## III The decay of global quantities

A first hint on the restoration of isotropy at large times can be obtained by the two-dimensional snapshots in Fig. 1. After a few eddy turnover times, it is evident that large-scale fluctuations become more and more isotropic. To give a quantitative measure, we collect for each run the temporal behavior of the following one-point quantities:

$`E_{ij}=\overline{v_i(t)v_j(t)},`$ (2)
$`\mathrm{\Omega }_{ij}=\overline{\omega _i(t)\omega _j(t)}.`$ (3)
By $`\overline{\mathrm{}}`$ we denote the average over space coordinates only, whereas $`\mathrm{}`$ will indicate the average over both initial conditions and space. The symmetric matrices $`E_{ij}(t)`$ and $`\mathrm{\Omega }_{ij}(t)`$ are then diagonalized at each time-step and the eigenvalues $`E_1(t),E_2(t),E_3(t)`$ and $`\mathrm{\Omega }_1(t),\mathrm{\Omega }_2(t),\mathrm{\Omega }_3(t)`$ are extracted. Since the forcing points in a fixed direction two eigenvalues are almost degenerate, say $`E_2`$ and $`E_3`$, and strongly differ from the first one, $`E_1`$. The typical decay of $`E_i(t)`$ and $`\mathrm{\Omega }_i(t)`$ for $`i=1,\mathrm{},3`$ is shown in Fig. 2. During the self-similar stage, $`t[10,10^6]`$, the energy eigenvalues fall off as $`E_{\{1,2,3\}}t^2`$, as expected for the decay in a bounded domain . The enstrophy eigenvalues, $`\mathrm{\Omega }_{\{1,2,3\}}`$ decay as $`t^{12/5}`$. The dimensional argument that captures these algebraic laws proceeds as follows. The energy decay is obtained by the energy balance $`dE(t)/dt=ϵ(t)`$, where we estimate $`Ev_{rms}^2(t)`$ and $`ϵ(t)v_{rms}^3(t)/L_0E^{3/2}(t)/L_0`$, and obtain $`E(t)t^2`$ and $`ϵ(t)t^3`$. As for the vorticity decay, we have $`\mathrm{\Omega }(t)(\delta _\eta v)^2/\eta ^2`$, where $`\eta `$ is the dissipative lengthscale and $`\delta _\eta v`$ is the typical velocity difference at separation $`\eta `$. Assuming a Kolmogorov scaling $`\delta _\eta vϵ(t)^{1/3}\eta ^{1/3}`$, and recalling that for a second-order hyperviscous dissipation $`ϵ(t)\nu (\delta _\eta v)^2/\eta ^4`$ we obtain $`\eta t^{3/10}`$ and $`\delta _\eta vt^{9/10}`$, whence $`\mathrm{\Omega }(t)t^{12/5}`$. We now focus on the process of recovery of isotropy in terms of global quantities. We identify two set of observables

$`\mathrm{\Delta }_{ij}E(t)`$ $`=`$ $`{\displaystyle \frac{E_i(t)E_j(t)}{E_i(t)+E_j(t)}},`$ (4)
$`\mathrm{\Delta }_{ij}\mathrm{\Omega }(t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_i(t)\mathrm{\Omega }_j(t)}{\mathrm{\Omega }_i(t)+\mathrm{\Omega }_j(t)}},`$ (5)
which vanish for isotropic statistics. Their rate of decay is therefore a direct measurement of the return to isotropy. The energy matrix $`E_{ij}`$ is particularly sensitive to the large scales while small-scale fluctuations are sampled by $`\mathrm{\Omega }_{ij}`$. As it can be seen from Fig. 3, both large and small scales begin to isotropize after roughly one eddy turnover time and become fully isotropic (within statistical fluctuations) after $`100`$ eddy turnover times. However, small scales show an overall degree of anisotropy much smaller than the large scales.

## IV The decay of small-scale fluctuations

### A The self-preservation hypothesis for anisotropic decay

The observables which characterize the decay of small-scale velocity fluctuations are the longitudinal structure functions

$$S^{(n)}(𝐫,t)=\left[\left(𝒗(𝒙+𝐫,t)𝒗(𝒙,t)\right)\widehat{𝐫}\right]^n.$$
(6)
Those quantities depend not only on the modulus of the separation $`𝐫`$ but also on its orientation $`\widehat{𝐫}`$, since the velocity statistics is anisotropic. A method to systematically disentangle isotropic from anisotropic contributions in the structure functions is based on the irreducible representations of the SO(3) group . In this approach, the observables (6) are expanded on the complete basis of the eigenfunctions of the rotation operator. The SO(3) decomposition of scalar objects, such as structure functions, is obtained by projection on the spherical harmonics $`Y_{jm}(\widehat{𝒓})`$:

$$S^{(n)}(𝐫,t)=\underset{j=0}{\overset{\mathrm{}}{}}\underset{m=j}{\overset{j}{}}S_{jm}^{(n)}(r,t)Y_{jm}(\widehat{𝒓}).$$
(7)
Here, $`S_{jm}^{(n)}(r,t)`$ denotes the projection of the $`n`$-th order structure function on the $`(j,m)`$ SO(3) sector, with $`j`$ and $`m`$ labeling the total angular momentum and its projection in the direction $`\widehat{𝒛}`$, respectively. Another equivalent possibility is to look at the SO(3) decomposition of the PDF of the longitudinal velocity differences. In this case, denoting by $`𝒫(\mathrm{\Delta },𝐫;t)`$ the probability that the longitudinal incremement $`\delta _rv\left(𝒗(𝐫,t)𝒗(0,t)\right)\widehat{𝒓}`$ be equal to $`\mathrm{\Delta }`$, we may project $`𝒫(\mathrm{\Delta },𝐫;t)`$ on the SO(3) basis:

$$𝒫(\mathrm{\Delta },𝐫;t)=\underset{j=0}{\overset{\mathrm{}}{}}\underset{m=j}{\overset{j}{}}𝒫_{jm}(r,\mathrm{\Delta };t)Y_{jm}(\widehat{𝒓}).$$
(8)
The projection $`𝒫_{jm}(r,\mathrm{\Delta };t)`$ plays the role of an “effective PDF” for each single SO(3) sector. (It should be however remarked that only the isotropic probability density $`𝒫_{00}(r,\mathrm{\Delta };t)`$ has the property of being everywhere positive and normalized to unity with respect to the weight $`r^2/(4\pi )`$.) Indeed, the projections of the longitudinal structure function on any sector $`(j,m)`$ can be reconstructed from the corresponding $`𝒫_{jm}(r,\mathrm{\Delta };t)`$ by averaging over all possible $`\mathrm{\Delta }`$’s:

$$S_{jm}^{(n)}(r,t)=𝑑\mathrm{\Delta }\mathrm{\Delta }^n𝒫_{jm}(r,\mathrm{\Delta };t).$$
(9)
That establishes the equivalence between the decompositions (7) and (8).
The main points broached here are about the long-time properties of the SO(3) projections, $`S_{jm}^{(n)}(r,t)`$. A simple isotropic generalization of the self-preservation hypothesis (see, e.g. Ref. ) amounts to writing:

$$S_{jm}^{(n)}(r,t)=V_{jm}^{(n)}\left(t\right)f_{jm}^{(n)}\left(r/L_{jm}(t)\right).$$
(10)
Note that $`V_{jm}^{(n)}\left(t\right)`$ takes explicitly into account the fact that large-scale properties may depend in a nontrivial way on both $`(j,m)`$ and the order $`n`$. Furthermore, $`L_{jm}(t)`$ accounts for the possibility that the characteristic length scale depend on the SO(3) sector.
In analogy with the observations made in the stationary case  we postulate a scaling behavior

$$S_{jm}^{(n)}(r,t)a_{jm}^{(n)}(t)\left(\frac{r}{L_{jm}(t)}\right)^{\zeta _j^{(n)}}.$$
(11)
The time behavior is encoded in both the decay of the overall intensity, accounted by the prefactors $`a_{jm}^{(n)}(t)`$, and the variation of the integral scales $`L_{jm}(t)`$. The representation (11) is the simplest one fitting the initial time statistics for $`t=0`$ and agreeing with the evolution given by the self preservation hypothesis in the isotropic case. The power law behavior for $`f_{jm}^{(n)}\left(r/L_{jm}(t)\right)`$ can be expected only in a time-dependent inertial range of scales $`\eta (t)rL(t)`$. As for the exponents appearing in (11), their values are expectedly the same as in the stationary case. In the latter situation it has been shown that they are organized hierarchically according to their angular sector $`j`$  :

$$\zeta _0^{(n)}\zeta _1^{(n)}\mathrm{}\zeta _j^{(n)}\mathrm{}.$$
(12)
Since the isotropic sector has the smallest exponent, at any given time and for a given intensity of the fluctuation (selected by the value of $`n`$) we have a recovery of isotropy going to smaller and smaller scales. Yet, deviations of the scaling exponents from their dimensional expectations make the recovery at small-scales much slower than what predicted by dimensional analysis . Moreover, there are quantities which should vanish in an isotropic field and actually blow up as the scale decreases .
Concerning the time evolution, it seems difficult to disentangle the dependence due to the decay of $`a_{jm}^{(n)}(t)`$ from the one due to the growth of the integral scale $`L_{jm}(t)`$. Here, we note only that the existence of a running reference scale, $`L_{jm}(t)`$ introduces some non-trivial relations between the spatial anomalous scaling and the decaying time properties, and those relations might be subject to experimental verification. In our case, the fact that the initial condition has a characteristic lengthscale comparable with the box size, simplifies the matter. Indeed we expect that $`L_{jm}(t)L_0`$, and the decay is due only to the fall off of the global intensity $`a_{jm}^{(n)}(t)`$. Unfortunately, a shortcoming is that the width of the inertial range $`L_0/\eta (t)`$ shrinks monotonically in time, thereby limiting the possibility of precise quantitative statements.

### B Numerical results

An overall view of the SO(3) projections at all resolved scales and for all measured decay times is presented in Fig. 4 for $`n=2`$ and the isotropic, $`(j=0,m=0)`$, sector. The same quantities are presented in Fig. 5 for the most intense anisotropic sector $`(j=4,m=0)`$.

We notice that as time elapses the dissipative range erodes the inertial one, as a consequence of the growth of the Kolmogorov scale $`\eta (t)`$. We notice in passing that the two curves corresponding to $`t=0`$ and $`t=\tau _0`$ almost coincide, i.e. even small scales are unchanged despite their typical eddy turnover times is much smaller than $`\tau _0`$. This finding has some consequences that will be discussed at length in section V. A similar qualitative trend is displayed by the most intense anisotropic sector, $`(j=4,m=0)`$ shown in Fig. 5, even though oscillations at small scales spoil significantly the scaling properties at small separations $`r`$.

In order to assess the relative importance of isotropic and anisotropic contributions, we plot in Fig. 6 the SO(3) projections of some $`(j,m)`$ sectors for the order $`n=2`$ and a fixed time, $`t=\tau _0`$. Fig. 7 shows the same quantities as in Fig. 6 but at a later time, $`t=100\tau _0`$. Although the small-scale behavior of anisotropic sectors readily becomes rather noisy, the various contributions are organized hierarchically and the isotropic contribution is dominant, as expected.
Let us now analyze quantitatively the time-decay of the structure functions at a fixed separation. In Fig. 8 we show the long-time decay of the second and fourth-order moments on the isotropic and an anisotropic sector at $`r=80`$, within the inertial range.

We observe that the anisotropic sectors decay faster than the isotropic one, i.e. at a fixed scale there is a tendency towards the recovery of isotropy at large times. The relative rate of decay can be quantified by the observable

$$\mathrm{\Pi }_{jm}^{(n)}(r,t)\frac{S_{jm}^{(n)}(r,t)}{S_{0,0}^{(n)}(r,t)}t^{\chi _j^{(n)}}.$$
(13)
In Fig. 9 we show $`\mathrm{\Pi }_{jm}^{(n)}(r,t)`$ at $`r=80`$ for structure functions of order $`n=2,4,6`$ and for $`(4,0)`$, one of the most intense anisotropic sectors. In the inset we also plot the same quantities at the smaller scale $`r=40`$. All anisotropic sectors, for all measured structure functions, decay faster than the isotropic one. The measured slope in the decay is about $`\chi _j^{(n)}0.3`$ almost independent, within the statistical errors, of the order $`n`$. Note that these results agree with the simple picture that the time-dependence in (11) is entirely carried by the prefactors $`a_{jm}^{(n)}(t)`$ and the value of the integral scales $`L_{jm}(t)`$ is saturated at the size of the box. Indeed, by assuming that large-scale fluctuations are almost Gaussian we have that the leading time-dependency of $`a_{jm}^{(2n)}`$ is given by $`a_{jm}^{(2)}a_{00}^{(2n2)}`$. For the isotropic sector, $`a_{00}^{(2n)}(a_{00}^{(2)})^n`$, and plugging that in (13), we get: $`\mathrm{\Pi }_{jm}^{(n)}(r,t)a_{jm}^{(2)}(t)/a_{00}^{(2)}(t)t^\chi ^{}`$ with $`\chi ^{}0.3(\pm 0.1)`$ independent of $`n`$. The quality of our data is insufficient to detect possible residual effects due to $`L_{jm}(t)`$, which would make $`\chi _j^{(n)}`$ depend on $`n`$ and $`j`$ because of spatial intermittency.
The interesting fact that we measure decay properties of the anisotropic sectors which are almost independent of the order of the structure functions indicate that we must expect some non-trivial time dependence in the shape of the PDF’s $`𝒫_{jm}(r,\mathrm{\Delta };t)`$ for $`j>0`$. The most accurate way to probe the rescaling properties of $`𝒫_{jm}(r,\mathrm{\Delta };t)`$ in time is to compute the generalized flatness:

$$K_{jm}^{(n)}(r,t)\frac{S_{jm}^{(n)}(r,t)}{\left(S_{jm}^{(2)}(r,t)\right)^{\frac{n}{2}}}t^{\alpha _j^{(n)}}$$
(14)
Were the PDF projection in the $`(j,m)`$ sector self-similar for $`t\tau _0`$, then $`K_{jm}^{(n)}(r,t)`$ would tend to costant values. This is not the case for anisotropic fluctuations, as it is shown in Fig. (10). The curves $`K_{jm}^{(n)}(r,t)`$ are collected for two fixed inertial range separations, $`r=80`$ and $`r=40`$ (inset), for two different orders, $`n=4,6`$ and for both the isotropic and one of the most intense anisotropic sectors$`(j=4,m=0)`$ . The isotropic flatness tends toward a constant value for large $`t`$. Conversely, its anisotropic counterparts are monotonically increasing with $`t`$, indicating a tendency for the anisotropic fluctuations to become more and more intermittent as time elapses. A consequence of the monotonic increase of intermittency for large times is the impossibility to find a rescaling function, $`g(t,r)`$, which makes the rescaled PDF $`g(t,r)𝒫_{jm}(r,\frac{\mathrm{\Delta }}{g(t,r)};t)`$ time-independent at large times. Let us notice that also the behavior in Fig. 10 is in qualitative agreement with the observation previously made that all time dependencies can be accounted by the prefactors $`a_{jm}^{(n)}(t)`$. Indeed, assuming that the length scales $`L_{jm}(t)`$ have saturated and that the large scale PDF is close to Gaussian, it is easy to work out the prediction $`K_{jm}^{(n)}(r,t)t^{\chi ^{}(1n/2)}`$, i.e. $`\alpha _j^{(n)}=\chi ^{}(n/21)`$.

We conclude this section by a brief summary of the results. We have found that isotropic fluctuations persist longer than anisotropic ones, i.e. there is a time-recovery, albeit slower than predicted by dimensional arguments, of isotropy during the decay process. We have also found that isotropic fluctuations decay in an almost self-similar way whilst the anisotropic ones become more and more intermittent. Qualitatively, velocity configurations get more isotropic but anisotropic fluctuations become, in relative terms, more “spiky” than the isotropic ones as time elapses.

## V Short-time decay

Let us now move to the properties of decay at short times ($`t\tau _0`$). Universality of small-scale forced turbulence is at the forefront of both theoretical and experimental investigation of real turbulent flows . The problem is to identify those statistical properties which are robust against changes of the large-scale physics, that is against changes in the boundary conditions and the forcing mechanisms. Our goal here is to relate the small-scale universal properties of forced turbulent statistics to those of short-time decay for an ensemble of initial configurations. An immediate remark is that one cannot expect an universal behaviour for all statistical observables as the very existence of anomalous scaling is the signature of the memory of the boundaries and/or the external forcing throughout all the scales. Indeed, the main message we want to convey here is that only the scaling of both isotropic and anisotropic small-scale fluctuations is universal, at least for forcings concentrated at large scales. The prefactors are not expected to be so. There is therefore no reason to expect that quantities such as the skewness, the kurtosis and in fact the whole PDF of velocity increments or gradients be universal.
This is the same behavior as for the passive transport of scalar and vector fields (see  and references therein). For those systems both the existence and the origin of the observed anomalous scaling laws have been understood and even calculated analytically for some instances in the special class of Kraichnan flows . Here, it is worth stressing that the universal character of scaling exponents is shared by both isotropic and anisotropic fluctuations .
For the Navier-Stokes case we have a huge amount of experimental and numerical indications that the velocity field shows anomalous scaling , suggesting the existence of phenomena “similar” to those of the linear case. However, carrying over the analytical knowledge developed for linear hydrodynamical problems involve some nontrivial, yet missing, steps. For the Navier-Stokes dynamics, linear equations of motion surface again but at the functional level of the whole set of correlation functions. In a schematic form:

$$_tC^{(n)}=\mathrm{\Gamma }^{(n+1)}C^{(n+1)}+\nu D^{(n)}C^{(n)}+F^{(n)},$$
(15)
where $`\mathrm{\Gamma }^{(n+1)}`$ is the integro-differential linear operator coming from the inertial and pressure terms and $`C^{(n+1)}`$ is a shorthand notation for a generic $`(n+1)`$-point correlator. The molecular viscosity is denoted by $`\nu `$ and $`D^{(n)}`$ is the linear operator describing dissipative effects. Finally, $`F^{(n)}`$ is the correlator involving increments of the large-scale forcing $`𝒇`$ and of the velocity field. Balancing inertial and injection terms gives dimensional scaling, and anomalously scaling terms must therefore have a different source. A natural possibility is that a mechanism similar to the one identified in linear transport problems be at work in the Navier-Stokes case as well. The anomalous contributions to the correlators would then be associated to statistically stationary solutions of the unforced equations (15). The scaling exponents would a fortiori be independent of the forcing and thus universal. As for the prefactors, the anomalous scaling exponents are positive and thus the anomalous contributions grow at infinity. They should then be matched at the large scales with the contributions coming from the forcing to ensure that the resulting combination vanish at infinity, as required for correlation functions. Our aim here is not to prove the previous points but rather to check over the most obvious catch: the Navier-Stokes equations being integro-differential, non-local contributions might directly couple inertial and injection scales and spoil the argument. This effect might be particularly relevant for anisotropic fluctuations where infrared divergences may appear in the pressure integrals .

In order to investigate the previous point, we performed two sets of numerical experiments in decay. The first set, A, is of the same kind as in the previous section, i.e. we integrated the unforced Navier-Stokes equations (1) with initial conditions picked from an ensemble obtained from a forced anisotropic stationary run. Statistical observables are measured as an ensemble average over the different initial conditions, $`_{ens}`$. The ensemble at the initial time of the decay process is therefore coinciding with that at the stationary state in forced runs. If correlation functions are indeed dominated at small scales by statistically stationary solutions of the unforced equations then the field should not decay. Specifically, the field should not vary for times smaller than the large-scale eddy turnover time $`\tau _0L_0/v^2^{1/2}`$, with $`L_0`$ denoting the integral scale of the flow. Those are the times when the effects of the forcing terms start to be felt. Note that this should hold at all scales, including the small ones whose turnover times are much faster than $`\tau _0`$. The second set of numerical simulations (set B) is meant to provide for a stringent test of comparison. The initial conditions are the same as before but for the random scrambling of the phases : $`\widehat{v}_i(𝒌)P_{il}(𝒌)\widehat{v}_l(𝒌)\mathrm{exp}(i\theta _l(𝒌))`$. Here, $`\widehat{v}`$ denotes the Fourier transform and $`P_{il}(𝒌)`$ is the incompressibility projector. In this way, the spectrum and its scaling are preserved but the wrong organization of the phases is expected to spoil the statistical stationarity of the initial ensemble. As a consequence, two different decays are expected for the two sets of experiments. In particular, contrary to set A, set B should vary at small scales on times of the order of the eddy turnover times $`\tau _rr^{2/3}`$. This is exactly what we found in the numerical simulations, as can be seen in Fig. 11, where the temporal behavior of longitudinal structure functions of order 2 and 4 is shown. The scaling of the contributions responsible for the observed behaviour at small scales are thus forcing independent.
As for anisotropic fluctuations, we also found two very different behaviors depending on the set of initial conditions. In Fig. 12 we show the case of the projection $`S_{jm}^{(n)}(r=60,t)`$ for the anisotropic sector $`j=4`$,$`m=0`$. As it can be seen, for set A of initial conditions the function is indeed not decaying up to a time of the order of $`\tau _0`$.
To conclude, the data presented here support the conclusion that nonlocal effects peculiar to the Navier-Stokes dynamics do not spoil arguments on universality based on analogies with passive turbulent transport. The picture of the anomalous contributions to the correlation functions having universal scaling exponents and non-universal prefactors follows.

## VI Conclusions

We have presented a numerical investigation of the decay of three-dimensional turbulence, both isotropic and anisotropic. Concerning short-time decay, we have compared the decay for two different sets of initial conditions, with and without phase correlations. That gave some new hints on the properties of universality of isotropic and anisotropic forced turbulence. As for long times, we have found that fluctuations in the inertial range become more and more isotropic. On the other hand, the anisotropic components become more and more intermittent, i.e. relatively intense anisotropic fluctuations become more and more probable. The main issue here was to investigate the effects of a bounded domain, i.e. situations where the integral scale cannot grow indefinitely. The decay process at long times is then governed by the set-up at large scales. Anisotropies decay in time at a rate almost independent of the order of the moments and of the kind of anisotropic fluctuations. Projections of the PDF’s on different SO(3) sectors show different intermittent properties (with the anisotropic sectors being more intermittent). An obvious further development of this study would be to investigate the case where the integral scales $`L_{jm}(t)`$ vary in time. Additional intermittency in time might then be brought in by the anomalous scaling in the space variables of the correlation functions.

This research was supported by the EU under the Grants No. HPRN-CT 2000-00162 “Non Ideal Turbulence” and No. HPRN-CT-2002-00300 ”Fluid Mechanical Stirring and Mixing”, and by the INFM (Iniziativa di Calcolo Parallelo). A.L. acknowledges the hospitality of the Observatoire de la Côte d’Azur, where part of the work has been done.


format: text

text:
# Stable 1:2 Resonant Periodic Orbits in Elliptic Three-Body Systems

## 1 Introduction

It is well known that orbital or mean-motion resonances are of great importance in the dynamics of planetary and satellite systems. Resonances drive the dynamical evolution of such systems by creating regions of stability and instability in orbital phase space. Mean-motion resonances are usually associated with orbits whose geometrical configurations are periodically repeated. Such orbits, known as resonant periodic orbits (hereafter RPOs), are of particular interest since they define the structure of the associated resonance. The stability of an RPO implies that it may be a potential location for harboring an actual body. Unstable RPOs, on the other hand, provide pathways for dramatic orbital evolution of bodies and can also be used to explain chaotic motions (Varadi $`\&`$ Kaula, 1999).

Until a few years ago, mean-motion resonances were found only in our solar system. Recent success of the precise radial velocity searches in detecting more than 100 extrasolar planets has extended the applicability of resonances to a broader context, several parsecs beyond the boundaries of our solar system.111See exoplanets.org for a complete and up-to-date list of extarsolar planets with their corresponding references. Among these planetary systems, GJ876 is a three-body system with two planets in a 1:2 resonance (Marcy et al., 2001) and 47 UMa has two planets that are a near 2:5 commensurability (Rivera $`\&`$ Haghighipour, 2003). Because of the universal importance and applicability of resonances, it is of great value to study RPOs and investigate their dynamical stability.

The most widely known RPOs are the ones associated with the 1:1 resonance. In a three-body system, these RPOs are the Lagrangian points in a frame co-rotating with the major planet. In case of the three-body system of Sun, Jupiter and a test particle, the stable RPOs are located at $`\pm 60^{}`$ with respect to the location of Jupiter on its orbit around the Sun. The roughly 1500 known Jupiter Trojan asteroids are the actual bodies that are trapped in the 1:1 resonance with Jupiter at these two locations.

Our knowledge of the RPOs of other orbital resonances is limited to only a few values of the eccentricity of the inner planet and its mass relative to the star. It is therefore essential to carry out studies similar to those for the 1:1 resonance for other mean-motion resonances and to explore their phase-parameter spaces for regions of stability. Such a study has been previously carried out by Varadi (1999) for the case of the 2:3 exterior resonance. In that paper, the method of differential continuation was employed to identify the regions of the phase-parameter space of planar, elliptic restricted three-body systems that correspond to stable orbits of the outer body while in the 2:3 resonance.

In this paper, we present the results of an extensive study of the RPOs of the exterior 1:2 resonance. The system of interest is planar, elliptic and restricted, with an inner massive planet (hereafter, the major planet) and an outer zero-mass particle. We employ the methodology of Varadi (1999) and vary the parameters of the system in a systematic way in order to identify its stable RPOs and their corresponding regions of stability in phase-parameter space.

A review of our methodology is presented in $`\mathrm{\S }`$ 2. Section 3 presents the results of our numerical study and in $`\mathrm{\S }`$ 4 we conclude by reviewing the results and discussing their implications to extrasolar planetary systems.

## 2 Methodology

Our method of identifying an RPO consists of two steps; search and continuation. To describe this process, we start by briefly reviewing the condition for an orbit to be periodic. Let us consider a planar, restricted three-body system with a mass-less particle in an exterior resonance with a major planet. If $`𝐱`$ denotes the four-dimensional vector of the initial position and velocity of the outer particle, then after a certain number of orbits of the major planet, the particle’s position and velocity will be given by another vector $`𝐅(𝐱)`$. The orbit is periodic if

$$𝐅(𝐱(\epsilon ),\epsilon )𝐱(\epsilon )=\mathrm{\hspace{0.17em}0},$$
(1)
where $`\epsilon `$ is a parameter of the system. For the systems considered here, $`\epsilon `$ can be taken to be the ratio of the mass of the major planet to that of the central star. The mapping $`𝐅`$ is given by numerically integrating the equations of motions. Solution vectors $`𝐱`$ of equation (1) represent RPOs of the system.

Equation (1) can be differentiated with respect to $`\epsilon `$ to obtain

$$\frac{d𝐱}{d\epsilon }=\left(\frac{d𝐅}{d𝐱}𝐈\right)^1\frac{d𝐅}{d\epsilon },$$
(2)
where $`𝐈`$ is the identity matrix. Equation (2) involves the variational equations of the dynamics of the system. It leads to a complicated system of differential equations in which the independent variable is $`\epsilon `$ (Varadi, 1999).

Equation (1) can, in principle, be solved by Newton’s method. Because the latter converges only for some vectors $`𝐱`$, due to the nonlinear nature of equation (1), many initial vectors have to be tried. For a chosen value of $`\epsilon `$, we select random initial conditions, integrate the dynamical equations and if $`𝐅`$ happens to be small, we apply Newton’s method to search for a nearby solution of equation (1). We then apply differential continuation (Keller, 1977), i.e., solve equation (2), to follow the changes in the RPO due to varying $`\epsilon `$.

We start our search for 1:2 resonant periodic orbits by choosing the semimajor axis of the outer body near its resonant value with the inner planet. The semimajor axis of the inner planet is assumed to be equal to unity. Many possible configurations are considered in which the other orbital parameters of the outer body are chosen randomly. The inner planet is placed at either the peri- or the apocenter without any loss of generality. We integrate the system for two orbits of the inner planet and form the differences between the initial and the final position and velocity vectors of the outer body. When these differences are small, we apply Newton’s method to further reduce them. This process is repeated until these differences become so small that the orbit can be considered to be an RPO.

The ratio of the mass of the inner planet to the mass of the central star $`(\mu )`$ and also the orbital eccentricity of the inner planet $`(e)`$ are the parameters of our system. During the process of differential continuation, we vary $`\mu `$ and $`e`$ and solve equation (2) to obtain new RPOs along a path in $`(\mu ,e)`$ space. To identify the regions of the phase-parameter space of the system which correspond to stable RPOs, we change $`\mu `$ and $`e`$ systematically and criss-cross the $`(\mu ,e)`$ space with paths of continuation along which the stability of RPOs is computed.

The stability of RPOs is examined by computing the eigenvalues of the monodromy matrix $`d𝐅/d𝐱`$ (Danby, 1964; Golub $`\&`$ Van Loan, 1996). Stable RPOs are centers of libration for nearby orbits. Unstable RPOs are associated with separatrix-like features (Guckenheimer $`\&`$ Holmes, 1983). The eigenvalues of the monodromy matrix provide the frequency of libration associated with stable RPOs and also the Lyapunov exponent of the outer body when the RPO is unstable. The system is linearly (spectrally) stable if these eigenvalues are on the unit circle in the complex plane. There is usually a separation of time scales between fast and slow dynamics and both can correspond to either stable or unstable motions. The fast dynamics usually consists of motions which are also present in the case of a circular orbit for the major planet. When this fast motion is stable, we call it libration. The slow dynamics is associated with the effects of orbital eccentricity of the major planet. We call stable motions in this case secular librations and unstable ones secular instabilities. An RPO is unstable when either of the two motions is unstable.

## 3 Numerical Results

We studied the orbit of the outer body of our three-body system for different values of the mass-ratio $`\mu `$ and also different values of the orbital eccentricity of the inner planet, $`e`$. To describe the shape and the orientation of these orbits, we use

$`h=e\mathrm{cos}\varpi ,k=e\mathrm{sin}\varpi `$ (3)
$`h^{}=e^{}\mathrm{cos}\varpi ^{},k^{}=e^{}\mathrm{sin}\varpi ^{},`$ (4)
where $`e^{}`$ is the orbital eccentricity of the outer body and $`\varpi `$ and $`\varpi ^{}`$ represent the longitudes of the pericenters of the inner planet and the outer body, respectively. In order to completely specify an RPO, one needs the $`h^{}`$, $`k^{}`$, semi-major axis and mean longitude of the RPO at a given orbital phase of the inner planet.

We performed an extensive initial numerical search for 1:2 resonant periodic orbits for several mass-ratios and orbital eccentricities of the inner planet. Once an RPO was found, we used its orbital parameters for starting the continuation process. Figure 1 shows the overall results of the continuation process for starting values of 0.001 and 0.1 for $`\mu `$ and $`e`$, respectively. For all the shown branches of continuation, we have we have $`\varpi `$ = 0 and therefore $`k=0`$ and $`h=e`$. For the same orbit of the inner planet, there can be several RPOs with different orbital elements on different branches. These branches are labeled according to the initial values of their $`h^{}`$ and $`k^{}`$ quantities. This method of labeling enables one to readily calculate the orbital eccentricity and the longitude of the pericenter of the outer body at the start of the continuation process. For instance, the branch labeled as $`h^{}52k^{}16`$ represents the result of continuation started from a resonant periodic orbit with $`h^{}=0.52`$ and $`k^{}=0.16`$. Such an orbit has an eccentricity equal to 0.544 and a longitude of pericenter equal to 17.1.

One can see from Figure 1 that certain branches shown here appear to be connected, either through so-called turning points or via circular connection. To explain these cases, we recall that any continuation branch corresponds to a periodic orbit that has been continued from an initial RPO, for different values of the mass and the orbital eccentricity of the inner body. The continuation process can break down at points where the matrix $`(d𝐅/d𝐱𝐈)`$ becomes singular. When two branches of continuation meet at a point on the $`h^{}`$-$`h`$ plane where $`h0`$, they are said to be connected through a turning point. Turning point is a general term indicating that a branch changes its direction, or in other word, at a turning point, the tangent line to the continuation branch will be vertical (e.g., Keller 1977). The two branches $`h^{}15k^{}00`$ and $`h^{}59k^{}00`$ are connected via a turning point. Figure 2 shows the changes of the different orbital parameters of the outer body for these two branches.

On the line $`h=0`$, the matrix above is always singular due to the fact that the differential equations for the outer planet become time-independent. In this case, the variational equations transport the equations of motions invariantly along the orbit (e.g., Siegel $`\&`$ Moser, 1971; Varadi, 1999). A circular connection occurs when two branches of continuation meet on the $`h=0`$ line. Branches $`h^{}004k^{}00`$ and $`h^{}15k^{}00`$, $`h^{}59k^{}00`$ and $`h^{}`$67$`k^{}00`$, and also $`h^{}`$05$`k^{}`$57 and $`h^{}52k^{}16`$ are connected in this way. We also have to note that the matrix above can also become singular at bifurcation points, which is discussed in the context of stability analysis.

The continuation process, however, can sometimes pass through singularities without major difficulties. This is due to the fact that the particular singularities are handled gracefully by the numerical integrator used to solve equation (2). An example is the smooth connection between $`h^{}004k^{}00`$ and $`h^{}15k^{}00`$ in Figure 1. Eventually there is a singularity of some type through which continuation is not possible.

As mentioned above, periodic orbits of the outer body are labeled by their orbital elements at the pericenter passage of the inner planet. In the cases in which the continuation process can be extended beyond $`h=0`$, i.e., beyond the orbit of the inner planet being circular, the pericenter and the apocenter of the inner planet exchange places. That is, pericenter becomes apocenter and vice versa. To assure that the orbital parameters of the outer body will still be calculated at the time of the pericenter passage of the inner planet, we integrate the equations of motion of the outer body for half of the orbital period of the inner planet and then rotate the coordinate systems in such a way that the pericenter of the inner planet is located in the positive direction of the $`x`$-axis. At the same time, since we are dealing with the 1:2 external resonance, the outer body travels only one-fourth of its orbit and one obtains a new set of orbital elements. Once the re-alignment process is completed, there will be two sets of orbital parameters associated with the outer body. We assure that these two sets describe orbits unambiguously by limiting the values of the mean longitude of the outer body, $`\lambda ^{}`$. In this study $`90^{}<\lambda ^{}<90^{}`$. The re-alignment correction becomes necessary when the orbital parameters of the outer body undergo sudden changes during the passage of the inner planet through a circular orbit. Figure 3 shows such changes for the re-aligned and non-realigned values of $`h`$ and $`h^{}`$ for two circularly connected branches of $`h^{}004k^{}00`$ and $`h^{}15k^{}00`$.

As mentioned earlier, Figure 1 provides an overall view of the different continuation branches for a specific value of the inner planet’s mass. When the latter is also varied, the branches become leaves parameterized by the mass and the eccentricity of the inner planet. In order to explore these leaves, the $`(\mu ,e)`$ space is criss-crossed with continuation paths. Along these paths, $`\mu `$ and $`e`$ are assumed to depend on a common independent variable. The derivatives of $`\mu `$ and $`e`$ with respect to the independent variable are prescribed to yield the desired continuation path. Since the actual variable used in the continuation is not $`e`$ but the velocity of the inner planet at peri- or apocenter, the paths will have slight curvature (see also Varadi, 1999).

We also studied the phase-parameter space of the system, in search of regions corresponding to stable resonant periodic orbits, for all different branches of Figure 1. In most cases, the RPOs were unstable. Figure 4 shows a criss-crossing of the $`(\mu ,e)`$ space with continuation paths of the RPO designated as $`h^{}004k^{}00`$. An interesting feature of this case is that the RPO can be stable in the horizontal directions while unstable in the vertical one. We note that planar motions, such as the RPOs in this study, can be unstable to vertical perturbations. Straightforward algebra reveals that the variational equations decouple the linearized dynamics of horizontal and vertical motions. Hence it is relatively simple to analyze vertical stability.

In Figure 5 and others of the same type, the quantities on the vertical axes represent either frequencies of libration or Lyapunov exponents, both computed from the eigenvalues of the monodromy matrix. When the plotted value is positive, it is known as the libration frequency. When the plotted values is negative, its absolute value is conventionally called Lyapunov exponent. Since both quantities have units of time, it is possible to plot them on the same axis. Varadi (1999) provides some clarifying examples. The orbit is linearly stable if all plotted values are positive and it is unstable otherwise. Since the linearized equations are uncoupled for horizontal and vertical motions, there is a single frequency of libration or Lyapunov exponent associated with the vertical motions. The horizontal motions are associated with a pair of quantities, each being either frequency of libration or Lyapunov exponent.

In Figure 5, we have plotted frequencies of librations or Lyapunov exponents along the vertical continuation path $`e=0.5`$ in Figure 4. For the values of the mass-ratio, $`\mu `$, considered in Figure 4, the orbit is unstable. From Figure 5a, one can see that while this instability is sustained along the vertical direction and for all values of $`\mu `$, the orbit becomes horizontally stable when $`0.0091\mu 0.03`$ and $`0.099\mu 0.115`$.

Figure 6 shows similar criss-crossing analysis for the continuation branch of the RPO $`h^{}`$-05$`k^{}`$-57. The system, in this case, reveals a region of stability with a distinct boundary from its large region of instability. The frequencies of librations or Lyapunov exponents of the system for different values of the mass-ratio $`\mu `$ and the orbital eccentricity of the inner planet are shown in Figures 7, 8 and 9. As shown in Figure 7, the orbit of the outer planet is stable along the horizontal line of $`\mu =0.001`$ in Figure 6, for all values of the orbital eccentricity of the outer planet less than 0.57. Along a vertical line such as $`e=0.5`$ in Figure 6, the system is stable for all values of the mass-ratio less than 0.0014 (Figure 8). Figure 9 shows frequencies of librations or Lyapunov exponents of the orbit of the outer body for the continuation branch in which the orbital eccentricity and the mass-ratio of the inner planet increase in steps of 0.01 and 0.001, respectively. The continuation path corresponding to this figure has been denoted by X in Figure 6. As shown in Figure 9, the orbit of the outer planet is stable for the values of the mass-ratio $`\mu `$ less than 0.005 and it becomes unstable for its larger values.

Figure 6 also shows that at higher values of the eccentricity of the inner planet, the stability of the system requires lower mass ratios implying that the RPOs of a system with a more massive inner planet in a 1:2 exterior resonance are most stable when the orbit of the inner planet is closer to a circular orbit. The actual boundaries of stable regions can be different for different resonances and different leaves of the same resonance. For instance, the case of the 2:3 external resonance has a small wedge of instability (Varadi, 1999).

Figures 10a and b show the different RPOs of such a system for the mass-ratio of 0.001 and different values of the inner planet’s orbital eccentricity. It is worth noting that for these stable orbits, the pericenter directions of the inner and outer planets are not aligned. The existence of a similar asymmetry of stable RPOs has also been suggested in the context of the Galilean satellites (Greenberg, 1987).

## 4 Discussion

We presented the results of a systematic search for stable resonant periodic orbits of planar, elliptic and restricted three-body systems in the exterior 1:2 resonance. Our approach in this paper toward the understanding of orbital resonances provides fundamental information on their structures. Stable periodic orbits are centers of librations whose periods are also obtained, at least for small amplitudes, from our stability computations.

Using the method of differential continuation, we examined the phase-parameter space of the system for different values of the mass-ratio and the orbital eccentricity of the inner planet. In the majority of the cases, the RPOs were unstable. We were, however, able to show that stable resonant periodic orbits can be found when the mass and the orbital eccentricity of the inner planet are within a certain range (Figure 6).

Our initial search for RPOs also covered the internal 2:1 resonance. We have not found any RPOs, either stable or unstable, in this case. This can have significant implications for the 2:1 Kirkwood gap in the asteroid belt. It might be possible to arrive at an explanation of the gap in terms of the fundamental properties of the Sun-asteroid-Jupiter three-body system. This approach would be an alternative to the one adopted by Moon, Morbidelli $`\&`$ Milgiorini (1998) which employ elaborate physical models.

The apparent lack of RPOs in the internal 2:1 resonance also raises questions regarding the same resonance in the full three-body system, i.e., when the mass of the outer body is non-zero. As the latter increases and the mass of the inner body decreases, the system passes from external resonance to internal. Along the way, one should find the limits of masses for which stable RPOs exist. Such studies will have immediate applications to some of the recently discovered extrasolar planetary systems such as GJ876 in which two planets are locked in a near 1:2 resonance.

We benefited from stimulating conversations with S. Peale. This work was partially supported by NASA Astrobiology Institute under grants NCC2-1056 (N.H.) and NCC2-1050 (J.C. and F.V.) and also by NASA NAG grant 5-11691 (F.V. and W.B.M.).


text:
# ’t Hooft Determinant: Fluctuations and Multiple Vacua

## I Introduction

The global $`U_L(3)\times U_R(3)`$ chiral symmetry of the QCD Lagrangian (for massless quarks) is broken by the $`U_A(1)`$ Adler-Bell-Jackiw anomaly of the $`SU(3)`$ singlet axial current $`\overline{q}\gamma _\mu \gamma _5q`$. Through the study of instantons , it became clear that effective $`2N_f`$ quark interactions, known as ’t Hooft interactions, violate $`U_A(1)`$, but are still invariant under chiral $`SU(N_f)\times SU(N_f)`$.In the case of two flavors they are four-fermion interactions, and the resulting low-energy theory resembles the old Vaks Larkin Nambu Jona-Lasinio model . In the case of three flavors they are six-fermion interactions which are responsible for the correct description of $`\eta `$ and $`\eta ^{}`$ physics, and additionally lead to the OZI-violating effects ,

$$_{\text{det}}=\kappa (\text{det}\overline{q}P_Rq+\text{det}\overline{q}P_Lq)$$
(1)
where the matrices $`P_{R,L}=(1\pm \gamma _5)/2`$ are projectors and the determinant is over flavor indices.

The physical degrees of freedom of QCD at low-energies are mesons. The bosonization of the effective quark interaction (1) by the path integral approach has been considered in , where the lowest order stationary phase approximation (SPA) has been used to estimate the leading contribution from the ’t Hooft determinant. In this approximation the functional integral is dominated by the stationary trajectories $`r_{\text{st}}(x)`$, determined by the extremum condition $`\delta S(r)=0`$ of the action $`S(r)`$. The lowest order SPA corresponds to the case in which the integrals associated with $`\delta ^2S(r)`$, for the path $`r_{\text{st}}(x)`$ are neglected and only $`S(r_{\text{st}})`$ contributes to the generating functional. The subject of the present work is to include the contribution associated with $`\delta ^2S(r)`$. Our motivation is the following: let us consider the one-dimensional integrals of the form

$$F(\lambda )=_{\mathrm{}}^{\mathrm{}}𝑑t\mathrm{exp}[i\lambda f(t)]$$
(2)
What is sought by the method of SP is the dominant contribution to $`F`$ as $`\lambda \mathrm{}`$. The dominant contribution to the integral comes from regions of $`t`$ where $`f^{}`$ vanishes. Supposing that $`f^{}`$ vanishes at only one point $`t_0`$ and neglecting contributions to the integral from regions far from $`t_0`$, one can obtain the result

$$F(\lambda )=\sqrt{\frac{2\pi i}{\lambda f^{\prime \prime }(t_0)}}e^{i\lambda f(t_0)}.$$
(3)
Hence $`F(\lambda )`$ goes to zero like $`1/\sqrt{\lambda }`$. The lowest order SPA would correspond to the result $`F(\lambda )=const\mathrm{exp}[i\lambda f(t_0)]`$, which has the incorrect asymptotic behavior $`F(\lambda )=𝒪(1)`$. If results related to finite-dimensional integrals, such as $`F(\lambda )`$, mean anything with regard to corresponding functional integrals, one can conclude that the leading SPA should include the contribution from the integrals associated with the second functional derivative $`\delta ^2S(r_{\text{st}})`$.

## II Path Integral Bosonization

To be definite, let us consider the theory of the quark fields in four dimensional Minkowski space, with dynamics defined by the Lagrangian density

$$=_{\text{NJL}}+_{\text{det}}.$$
(4)
The first term here is the extended version of the Nambu – Jona-Lasinio (NJL) Lagrangian $`_{\text{NJL}}=_0+_{\text{int}}`$, consisting of the free field part

$$_0=\overline{q}(i\gamma ^\mu _\mu \widehat{m})q,$$
(5)
and the $`U(3)_L\times U(3)_R`$ chiral symmetric four-quark interaction

$$_{\text{int}}=\frac{G}{2}[(\overline{q}\lambda _aq)^2+(\overline{q}i\gamma _5\lambda _aq)^2].$$
(6)
We assume that quark fields have color and flavor indices running through the set $`i=1,2,3`$; $`\lambda _a`$ are the standard $`U(3)`$ Gell-Mann matrices with $`a=0,1,\mathrm{},8`$. The current quark mass, $`\widehat{m}`$, is a nondegenerate diagonal matrix with elements $`\text{diag}(\widehat{m}_u,\widehat{m}_d,\widehat{m}_s)`$, it explicitly breaks the global chiral $`U(3)_L\times U(3)_R`$ symmetry of the $`_{\text{NJL}}`$ Lagrangian. The second term in (4) is given by (1). Letting

$$s_a=\overline{q}\lambda _aq,p_a=\overline{q}i\gamma _5\lambda _aq,s=s_a\lambda _a,p=p_a\lambda _a$$
(7)
yields

$$_{\text{det}}=\frac{\kappa }{64}\left[\text{det}(s+ip)+\text{det}(sip)\right].$$
(8)
The dynamics of the system is described by the vacuum transition amplitude in the form of the path integral

$$Z=𝒟q𝒟\overline{q}\mathrm{exp}\left(id^4x\right).$$
(9)
which can be written in the equivalent form

$$Z=𝒟\sigma _a𝒟\varphi _a𝒟q𝒟\overline{q}\mathrm{exp}\left(id^4x_q\right)𝒟r_{1a}𝒟r_{2a}\mathrm{exp}\left(id^4x_r\right)$$
(10)
where

$`_q(\overline{q},q,\sigma ,\varphi )`$ $`=`$ $`\overline{q}(i\gamma ^\mu _\mu \widehat{m}\sigma +i\gamma _5\varphi )q,`$ (11)
$`_r(\sigma ,\varphi ,r_1,r_2)`$ $`=`$ $`2G\left[(r_{1a})^2+(r_{2a})^2\right]2(r_{1a}\sigma _a+r_{2a}\varphi _a)`$ (12)
$``$ $`{\displaystyle \frac{\kappa }{8}}\left[\text{det}(r_1+ir_2)+\text{det}(r_1ir_2)\right].`$ (13)
The Fermi fields enter the action bilinearly, we can always integrate over them. At this stage one should also shift the scalar fields $`\sigma _a\sigma _a+\mathrm{\Delta }_a`$ by demanding that the vacuum expectation values of the shifted fields vanish $`<0|\sigma _a|0>=0`$, yielding gap equations to fix parameters $`\mathrm{\Delta }_a=m_a\widehat{m}_a`$, where $`m_a`$ denotes the constituent quark masses ***The shift by the current quark mass is needed to hit the correct vacuum state, see e.g. .. To evaluate path integrals over $`r_{1,2}`$ one has to use the method of stationary phase, or, after the formal analytic continuation in the time coordinate $`x_4=ix_0`$, the method of steepest descents. The analytical continuation of the Euclidean version of the path integral under consideration is (see  for details),

$$J(\sigma ,\varphi )=_{i\mathrm{}+r_{\text{st}}}^{+i\mathrm{}+r_{\text{st}}}𝒟r_{1a}𝒟r_{2a}\mathrm{exp}\left(d^4x_r(\sigma ,\varphi ,r_1,r_2)\right).$$
(14)
Near the saddle point $`r_{\text{st}}^a`$,

$$_r_r(r_{\text{st}})+\frac{1}{2}\underset{\alpha ,\beta }{}\stackrel{~}{r}_\alpha _{\alpha \beta }^{\prime \prime }(r_{\text{st}})\stackrel{~}{r}_\beta $$
(15)
where $`r_{\text{st}}^a`$ is a solution of the equations $`_r^{}(r_1,r_2)=0`$ determining a flat spot of the surface $`_r(r_1,r_2)`$,

$$\{\begin{array}{ccc}\hfill 2Gr_1^a(\sigma +\mathrm{\Delta })_a\frac{3\kappa }{8}A_{abc}(r_1^br_1^cr_2^br_2^c)& =& 0\hfill \\ \hfill 2Gr_2^a\varphi _a+\frac{3\kappa }{4}A_{abc}r_1^br_2^c& =& 0.\hfill \end{array}$$
(16)
with totally symmetric constants, $`A_{abc}`$, closely related with the $`U(3)`$ constants $`d_{abc}`$. We use in (15) symbols $`\stackrel{~}{r}^a`$ for the differences $`(r^ar_{\text{st}}^a)`$. To deal with the multitude of integrals in (14) we define a column vector $`\stackrel{~}{r}`$ with eighteen components $`\stackrel{~}{r}_\alpha =(\stackrel{~}{r}_1^a,\stackrel{~}{r}_2^a)`$ and with the matrix $`_{\alpha \beta }^{\prime \prime }(r_{\text{st}})`$ being equal to

$$_{\alpha \beta }^{\prime \prime }(r_{\text{st}})=4GQ_{\alpha \beta },Q_{\alpha \beta }=\left(\begin{array}{cc}\delta _{ab}\frac{3\kappa }{8G}A_{abc}r_{1\text{st}}^c& \frac{3\kappa }{8G}A_{abc}r_{2\text{st}}^c\\ \frac{3\kappa }{8G}A_{abc}r_{2\text{st}}^c& \delta _{ab}+\frac{3\kappa }{8G}A_{abc}r_{1\text{st}}^c\end{array}\right).$$
(17)
Eq.(14) can now be concisely written as

$$J(\sigma ,\varphi )=\mathrm{exp}\left(d^4x_r(r_{\text{st}})\right)_i\mathrm{}^{+i\mathrm{}}𝒟\stackrel{~}{r}_\alpha \mathrm{exp}\left(2Gd^4x\stackrel{~}{r}^\text{t}Q(r_{\text{st}})\stackrel{~}{r}\right)\left[1+𝒪(\mathrm{})\right].$$
(18)
The solutions of Eqs.(16) are the following even and odd parity combinations $`r_{1\text{st}}^a`$ and $`r_{2\text{st}}^a`$ expressed in the form of increasing powers in $`\sigma _a,\varphi _a`$

$`r_{1\text{st}}^a`$ $`=`$ $`h_a+h_{ab}^{(1)}\sigma _b+h_{abc}^{(1)}\sigma _b\sigma _c+h_{abc}^{(2)}\varphi _b\varphi _c+\mathrm{}`$ (19)
$`r_{2\text{st}}^a`$ $`=`$ $`h_{ab}^{(2)}\varphi _b+h_{abc}^{(3)}\varphi _b\sigma _c+\mathrm{}`$ (20)
Putting these expansions in Eqs.(16) one obtains a series of selfconsistent equations to determine the constants $`h_a`$, $`h_{ab}^{(1)}`$, $`h_{ab}^{(2)}`$, etc.  As a result we get

$$_r(r_{\text{st}})=2h_a\sigma _ah_{ab}^{(1)}\sigma _a\sigma _bh_{ab}^{(2)}\varphi _a\varphi _b+𝒪(\text{field}^3).$$
(21)
We now turn to the evaluation of the path integral in Eq.(18). In order to define the measure $`𝒟\stackrel{~}{r}_\alpha `$ more accurately we expand $`\stackrel{~}{r}_\alpha `$ in a Fourier series

$$\stackrel{~}{r}_\alpha (x)=\underset{n=1}{\overset{\mathrm{}}{}}c_{n,\alpha }\phi _n(x),$$
(22)
assuming that suitable boundary conditions are imposed. The set of the real functions $`\{\phi _n(x)\}`$ form an orthonormal and complete sequence, therefore

$$𝒟\stackrel{~}{r}_\alpha \mathrm{exp}\left(2Gd^4x\stackrel{~}{r}^\text{t}Q(r_{\text{st}})\stackrel{~}{r}\right)=\frac{C}{\sqrt{det(2G\lambda _{nm}^{\alpha \beta })}}.$$
(23)
with

$$2G\lambda _{nm}^{\alpha \beta }=\left(\begin{array}{cc}h_{ac}^{(1)1}& 0\\ 0& h_{ac}^{(2)1}\end{array}\right)_{\alpha \sigma }\left(\delta _{\sigma \beta }\delta _{nm}+d^4x\phi _n(x)F_{\sigma \beta }(x)\phi _m(x)\right)$$
(24)
and

$$F_{\sigma \beta }=\frac{3\kappa }{4}A_{eba}\left(\begin{array}{cc}h_{ce}^{(1)}(r_{1\text{st}}^ah_a)& h_{ce}^{(1)}r_{2\text{st}}^a\\ h_{ce}^{(2)}r_{2\text{st}}^a& h_{ce}^{(2)}(r_{1\text{st}}^ah_a)\end{array}\right)_{\sigma \beta }.$$
(25)
Only the matrix $`F_{\sigma \beta }`$ depends here on fields $`\sigma ,\varphi `$. By absorbing in $`C`$ the irrelevant field independent part of $`2G\lambda _{nm}^{\alpha \beta }`$, and expanding the logarithm in the representation $`det(1+F)=\mathrm{exp}\text{tr}\mathrm{ln}(1+F)`$, one can obtain finally for the complete action

$$S_r=d^4x\left\{_r(r_{\text{st}})+\frac{a}{2G^2}\underset{n=1}{\overset{\mathrm{}}{}}\frac{(1)^n}{n}\text{tr}\left[F_{\alpha \beta }^n(r_{\text{st}})\right]\right\}$$
(26)
proposing that the undetermined dimensionless constant $`a`$ will be fixed by confronting the model with experiment afterwards

## III The Ground State

Let us study the ground state of the model under consideration, then properties of the excitations will follow naturally. In Eq.(12) the field coefficients $`h_i`$ obey

$$2Gh_i\mathrm{\Delta }_i=\frac{\kappa }{8}t_{ijk}h_jh_k$$
(27)
with the totally symmetric coefficients $`t_{ijk}`$ equal to zero except for $`t_{uds}=1`$ and with order parameters $`\mathrm{\Delta }_i0(i=u,d,s)`$.

A tadpole graphs calculation gives for the gap equations the following result

$$2h_i+\frac{3a\kappa }{8G^2}\left(h_{ab}^{(2)}h_{ab}^{(1)}\right)A_{abc}h_{ci}^{(1)}=\frac{N_c}{2\pi ^2}m_iJ_0(m_i^2)$$
(28)
where the left hand side is the contribution from (26) and the right hand side is the contribution of the quark loop from (10), .

The second term on the left hand side of Eq.(28) is the correction resulting from the Gaussian integrals of the steepest descent method, comprising the effects of small fluctuations around the stationary path. If one puts for a moment $`a=0`$ in Eq.(28), and combines the result with Eqs.(27), one finds gap equations which are very similar to the ones obtained in . For this case one obtains values of $`(m_u,m_s)`$ which are uniquely related to values of $`(G,\kappa )`$.

The general case, which we have when $`a>0`$ in Eq.(28), yields in turn a region for $`m_u,m_s`$, where three values of couplings $`(G,\kappa )`$ are possible.

Conversely, one can study the solutions: $`m_u=m_u(G,\kappa ),m_s=m_s(G,\kappa )`$ at fixed values of input parameters: $`\mathrm{\Lambda },\widehat{m}_u=\widehat{m}_d,\widehat{m}_s`$. Again we find several extremal solutions. The minima will be identified in a future analysis of the effective potential. We refer to  for graphical displays.

Multivalued solutions of the gap equations have been obtained in a different context in .

## IV Concluding remarks

The purpose of this work has been twofold. Firstly we have developed the technique which is necessary to go beyound the lowest order SPA in the problem of the path integral bosonization of the ’t Hooft six quark interaction. This technique is rather general and can be readily used in other applications. Second, we have explored with considerable detail the implications of taking the quantum fluctuations in account in the description of the hadronic vacuum. We encountered several classes of solutions. These multiple vacua may have very interesting physical consequences and applications.

Acknowledgemnets: This work is supported by Fundação para a Ciência e a Tecnologia, POCTI/35304/FIS/2000 and NATO ”Outreach” Cooperation Program.


text:
# Tailorable acceptor C_{60-𝑛}B_𝑛 and donor C_{60-𝑚}N_𝑚 pairs for molecular electronics

## Abstract

Our first-principles calculations demonstrate that $`\mathrm{C}_{60n}\mathrm{B}_n`$ and $`\mathrm{C}_{60m}\mathrm{N}_m`$ molecules can be engineered as the acceptors and donors, respectively, which are needed for molecular electronics, by properly controlling the dopant number $`n`$ and $`m`$ in $`\mathrm{C}_{60}`$. As an example, we show that acceptor $`\mathrm{C}_{48}\mathrm{B}_{12}`$ and donor $`\mathrm{C}_{48}\mathrm{N}_{12}`$ are promising components for molecular rectifiers, carbon nanotube-based $`n`$-$`p`$-$`n`$ ($`p`$-$`n`$-$`p`$) transistors and $`p`$-$`n`$ junctions.

PACS (numbers): 33.15.-e, 85.65.+h, 31.15.Ar

Modern microelectronics and computation are advancing at an extremely fast rate because of remarkable circuit miniaturization. However, this trend will soon reach the scale of atoms or molecules. To continue toward faster and smaller computers, new schemes are required. Molecular electronics is one such approach.

One major problem in molecular electronics is connecting the functional molecules and assigning the observed electrical properties in an unambiguous way to the molecules in question. Fullerenes are large enough to be identified by transmission electron microscopy or scanning probe methods, are stable and easy to build into molecular circuits , and might be inserted into single-walled carbon nanotubes (SWNT). Hence, fullerenes should be ideal components for molecular electronics.

As in semiconductor electronics, acceptor/donor pairs are critical for use in molecular electronics, for example, molecular rectifiers , nanoscale $`p`$-$`n`$-$`p`$ transistors and $`p`$-$`n`$ junctions. For traditional silicon doping, group V atoms (for example, phosphorous) act as donors and group III ( for example, boron) are acceptors. Analogous acceptor/donor schemes are needed in molecular electronics. Fullerenes are unique because they can be doped in several different ways (for example, endohedral, substitutional and exohedral doping). This should provide a wide range of possible acceptor/donor schemes.

To design active molecular devices, components with controllable electronic properties are needed. For example, to obtain molecular rectification, the lowest unoccupied molecular orbital (LUMO) of the acceptor should lie at or slightly above the Fermi level of the electrode and above the highest occupied molecular orbital (HOMO) of the donor. Hence it is important to search for desired acceptor/donor pairs which satisfy the requirement. In this letter, we suggest a controlled approach to obtain such pairs from $`\mathrm{C}_{60}`$ molecules by using substitutional doping. Because the average carbon-carbon bond length in $`\mathrm{C}_{60}`$ is slightly larger than that in graphite, which can only be substitutionally doped by boron, and the force constants are somewhat weakened by the curvature of the $`\mathrm{C}_{60}`$ surface, both boron and nitrogen can substitute for one or more carbons in $`\mathrm{C}_{60}`$. Our first-principles calculations demonstrate that $`\mathrm{C}_{60n}\mathrm{B}_n`$ and $`\mathrm{C}_{60m}\mathrm{N}_m`$ molecules can be engineered as the acceptors and donors, respectively, which are desired for molecular electronics by properly controlling the dopant number $`n`$ and $`m`$. As an example, we present the electronic properties of the acceptor $`\mathrm{C}_{48}\mathrm{B}_{12}`$ and donor $`\mathrm{C}_{48}\mathrm{N}_{12}`$, and discuss their potential applications in molecular electronics.

First, we discuss the electronic properties of $`\mathrm{C}_{60n}\mathrm{B}_n`$, $`\mathrm{C}_{60n}\mathrm{N}_n`$ and $`\mathrm{C}_{60}`$ to show that $`\mathrm{C}_{60n}\mathrm{B}_n`$ and $`\mathrm{C}_{60n}\mathrm{N}_n`$ can act as controlled dopants. The number $`n`$ ranges from 1 to 12. Since there are many isomers of $`\mathrm{C}_{60n}\mathrm{X}_n`$ at fixed $`n`$, we only take the dopant assignment of Hultman et al.  as an example: each pentagon receives a maximum of one dopant X, and two X are separated from each other by two carbon atoms. The optimized geometry and total energy are calculated by using the Gaussian 98 program with the B3LYP hybrid density functional theory (DFT) and 6-31G(d) basis set. The calculated results are summarized in Table I. It is found that the binding energy of $`\mathrm{C}_{60n}\mathrm{X}_n`$ (X=B,N) decreases monotonically with increasing integer $`n`$. For $`\mathrm{C}_{60n}\mathrm{B}_n`$, the binding energy is 0.21 eV/atom to 0.38 eV/atom lower than that ($`E_b=6.98`$ eV/atom) of $`\mathrm{C}_{60}`$. For $`\mathrm{C}_{60n}\mathrm{N}_n`$, it is 0.10 eV/atom to 0.61 eV/atom lower than that of $`\mathrm{C}_{60}`$. The stabilities of $`\mathrm{C}_{60n}\mathrm{B}_n`$ and $`\mathrm{C}_{60n}\mathrm{N}_n`$ are comparable to but less than that of $`\mathrm{C}_{60}`$. In 1991, Smalley and coworkers successfully synthesized $`\mathrm{C}_{60n}\mathrm{B}_n`$ with $`1n6`$. In 1995, Hummelen et al. produced $`\mathrm{C}_{59}\mathrm{N}`$. Very recently, Hultman et al. have reported the existence of $`\mathrm{C}_{48}\mathrm{N}_{12}`$. These experiments show that certain stable $`\mathrm{C}_{60n}\mathrm{X}_n`$ (X=B, N) can be made. Our results show that other B- or N-doped $`\mathrm{C}_{60n}\mathrm{X}_n`$ structures for $`n12`$ have similar stability. Our calculated ionization potential ($`E_I=7.32`$ eV) and electron affinity ($`E_A=2.40`$ eV) for $`\mathrm{C}_{60}`$ agree well with experiments ($`\mathrm{E}_I=7.54\pm 0.01\mathrm{eV}`$  and $`\mathrm{E}_A=2.689\pm 0.008\mathrm{eV}`$ ), indicating the accuracy of our calculations for $`\mathrm{C}_{60n}\mathrm{X}_n`$. Among $`\mathrm{C}_{60}`$ and $`\mathrm{C}_{60n}\mathrm{X}_n`$ (X=B,N) at fixed $`n`$, $`\mathrm{C}_{60n}\mathrm{B}_n`$ has the highest electron affinity, while $`\mathrm{C}_{60n}\mathrm{N}_n`$ has the lowest ionization potential, indicating that $`\mathrm{C}_{60n}\mathrm{B}_n`$ and $`\mathrm{C}_{60n}\mathrm{N}_n`$ can serve as electron acceptor and donor, respectively.

Calculations of LUMO and HOMO energies are necessary to explore combinations of donors and acceptors suitable for molecular electronics. Table I shows the change of the LUMO/HOMO energies (and LUMO-HOMO gaps $`E_g=E_lE_h`$) as we vary the dopant number $`n`$ in the $`\mathrm{C}_{60n}\mathrm{X}_n`$ molecule. Hence, acceptor/donor pairs which satisfy the required conditions in molecular electronics can be obtained by controlling the dopant number $`n`$ and $`m`$ in $`\mathrm{C}_{60n}\mathrm{B}_n`$ and $`\mathrm{C}_{60m}\mathrm{N}_m`$. For example, based on the acceptor LUMO and donor HOMO shown in Table I, we may choose acceptor/donor pairs, such as C58B2/C58N2, C54B6/C54N6, C49B11/C51N9, or C48B12/C48N12 to build molecular rectifiers exhibiting similar rectification behaviors. In the following, we take C48B12/C48N12 as one example to show their applications in molecular electronics.

Table I: LUMO ($`E_l`$, in eV), HOMO ($`E_h`$, in eV), binding energy ($`E_b`$, in eV/atom), ionization potential ($`E_I`$, in eV), and electron affinity ($`E_A`$, in eV) calculated for $`\mathrm{C}_{60n}\mathrm{X}_n`$ (X=B,N) using B3LYP/6-31G(d).

| n | $`\mathrm{C}_{60n}\mathrm{B}_n`$ | | | | | $`\mathrm{C}_{60n}\mathrm{N}_n`$ | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | $`E_l`$ | $`E_h`$ | $`E_b`$ | $`E_I`$ | $`E_A`$ | $`E_l`$ | $`E_h`$ | $`E_b`$ | $`E_I`$ | $`E_A`$ |
| 1 | -4.26 | -5.57 | 6.77 | 6.77 | 3.07 | -3.63 | -4.72 | 6.88 | 5.93 | 2.15 |
| 2 | -4.64 | -5.25 | 6.75 | 6.35 | 3.53 | -3.71 | -4.64 | 6.82 | 5.85 | 2.27 |
| 3 | -4.21 | -5.49 | 6.74 | 6.37 | 3.32 | -3.72 | -4.69 | 6.78 | 5.92 | 2.29 |
| 4 | -4.43 | -5.60 | 6.73 | 6.79 | 3.32 | -3.91 | -4.56 | 6.73 | 5.83 | 2.19 |
| 5 | -4.59 | -5.54 | 6.71 | 6.85 | 3.56 | -3.74 | -4.56 | 6.69 | 5.95 | 2.20 |
| 6 | -4.55 | -5.49 | 6.70 | 7.01 | 2.98 | -3.74 | -4.61 | 6.64 | 5.85 | 2.32 |
| 7 | -4.49 | -5.54 | 6.67 | 6.75 | 3.23 | -3.58 | -4.59 | 6.60 | 5.91 | 2.39 |
| 8 | -4.62 | -5.34 | 6.66 | 6.44 | 3.50 | -3.60 | -4.41 | 6.55 | 5.65 | 2.36 |
| 9 | -4.65 | -5.43 | 6.64 | 6.61 | 3.34 | -3.49 | -4.45 | 6.51 | 5.77 | 1.90 |
| 10 | -4.57 | -5.44 | 6.63 | 6.57 | 3.45 | -3.31 | -4.28 | 6.46 | 5.52 | 1.89 |
| 11 | -4.45 | -5.61 | 6.62 | 6.87 | 3.30 | -3.36 | -4.35 | 6.41 | 5.71 | 2.29 |
| 12 | -4.24 | -5.58 | 6.60 | 6.73 | 3.08 | -2.61 | -4.38 | 6.37 | 5.66 | 1.49 |

The molecular geometry (size) and symmetry play important roles in molecular electronics . For example, an important factor in inducing rectification is some geometric asymmetry in the molecular junction. $`\mathrm{C}_{60}`$ is a truncated icosahedron with a perfect $`I_h`$ symmetry. Then, the changes in geometry (size) and symmetry of $`\mathrm{C}_{60n}\mathrm{X}_n`$, due to the dopant-induced effects, are needed to know. Fig.1 shows the optimized structure of $`\mathrm{C}_{48}\mathrm{B}_{12}`$. For comparison, similar calculations were done for $`\mathrm{C}_{60}`$ and $`\mathrm{C}_{48}\mathrm{N}_{12}`$. The equilibrium $`\mathrm{C}_{48}\mathrm{B}_{12}`$, similar to $`\mathrm{C}_{48}\mathrm{N}_{12}`$, has one dopant (boron) per pentagon and two dopants (boron) preferentially sit in a hexagon. In $`\mathrm{C}_{48}\mathrm{X}_{12}`$ (X=B, N), the doping-induced distortion from the perfect $`\mathrm{C}_{60}`$ sphere is not localized to the neighborhood of each dopant but extends throughout the whole cage. The symmetry of $`\mathrm{C}_{48}\mathrm{B}_{12}`$ is the $`C_i`$ point group, different from the $`S_6`$ symmetry for $`\mathrm{C}_{48}\mathrm{N}_{12}`$. This is due to their opposite electronic polarization, as discussed below. The molecular structure is examined by calculating the average distances (or radii $`\mathrm{R}`$) from each atom to the center of the molecule. $`\mathrm{C}_{48}\mathrm{B}_{12}`$ has an ellipsoidal structure, as in $`\mathrm{C}_{48}\mathrm{N}_{12}`$, with 10 unique sites labelled 1 to 10 in Fig.1 and 10 different radii ranging from 0.347 nm (site 2, 2’) to 0.387 nm (site 8, 8’). In contrast, the 10 unique radii for $`\mathrm{C}_{48}\mathrm{N}_{12}`$ range from 0.340 nm (site 7) to 0.362 nm (site 5). For $`\mathrm{C}_{60}`$, each carbon atom has an equal radius $`\mathrm{R}=0.355\mathrm{nm}`$, the same as that found in experiment . Furthermore, we find 15 unique bonds (specified in Fig.1) in $`\mathrm{C}_{48}\mathrm{B}_{12}`$: 6 boron-carbon bonds with lengths between 0.154 nm and 0.159 nm, and 9 carbon-carbon bonds with 0.139 nm to 0.150 nm lengths. In comparison, $`\mathrm{C}_{48}\mathrm{N}_{12}`$ has 6 nitrogen-carbon bonds with lengths ranging from 0.141 nm to 0.143 nm, and 9 carbon-carbon bonds with lengths from 0.139 nm to 0.145 nm. $`\mathrm{C}_{60}`$, however, has one kind of single C-C bond (0.145 nm) and one kind of double C=C bond (0.139 nm), which are in excellent agreement with experiments (0.14459 nm and 0.13997 nm ).

FIG. 1.: $`\mathrm{C}_{48}\mathrm{B}_{12}`$ structure optimized with B3LYP/6-31G(d). Red (grey) is for B (C) sites. The 15 bonds between labeled vertices are all unique.

To determine the electronic polarization, Mulliken charge analysis  was performed for the on-site charge $`\mathrm{Q}_m`$ of atoms in $`\mathrm{C}_{48}\mathrm{B}_{12}`$, $`\mathrm{C}_{48}\mathrm{N}_{12}`$ and $`\mathrm{C}_{60}`$. $`\mathrm{C}_{48}\mathrm{B}_{12}`$ has two types of boron dopants: one boron with net Mulliken charge $`\mathrm{Q}_\mathrm{m}=0.164|e|`$ (site 5) and one with $`\mathrm{Q}_\mathrm{m}=0.187|e|`$ (site 9), while carbon atoms in $`\mathrm{C}_{48}\mathrm{B}_{12}`$ have negative $`\mathrm{Q}_\mathrm{m}`$ in the range from -0.004 $`|e|`$ (site 3) to -0.083 $`|e|`$ (site 8). In $`\mathrm{C}_{48}\mathrm{N}_{12}`$, there are two types of nitrogen dopants: one nitrogen with $`\mathrm{Q}_\mathrm{m}=0.595|e|`$ (site 5) and one with $`\mathrm{Q}_\mathrm{m}=0.600|e|`$ (site 9), and two types of carbon atoms: one-fourth of the carbon atoms with negative $`\mathrm{Q}_\mathrm{m}`$ (-0.0125 $`|e|`$ at site 2 and -0.0298 $`|e|`$ at site 3) and three-fourths of the carbon atoms with $`\mathrm{Q}_\mathrm{m}`$ in the range of 0.192 $`|e|`$ (site 6) to 0.227 $`|e|`$ (site 4). In $`\mathrm{C}_{60}`$, $`\mathrm{Q}_\mathrm{m}=0`$ for each carbon atom. Although the Mulliken analysis cannot estimate the atomic charges quantitatively, their signs are correct. Hence, $`\mathrm{C}_{60}`$ is isotropic, while $`\mathrm{C}_{48}\mathrm{B}_{12}`$ and $`\mathrm{C}_{48}\mathrm{N}_{12}`$ have opposite electronic polarizations, leading to their different symmetries.

FIG. 2.: Orbital energies $`E_k`$ and symmetries of the $`k`$th eigenstate for $`\mathrm{C}_{48}\mathrm{X}_{12}`$ (X=B,C,N) calculated with B3LYP/6-31G(d).

Applications of acceptors $`\mathrm{C}_{60n}\mathrm{B}_n`$ and donors $`\mathrm{C}_{60n}\mathrm{N}_n`$ require the proper line-up of energy levels. Upon boron or nitrogen doping, the degeneracy of the original $`\mathrm{C}_{60}`$ levels, shown in Fig.2, is removed by the structural distortion of the icosahedral symmetry. $`\mathrm{C}_{48}\mathrm{B}_{12}`$ ($`\mathrm{C}_{48}\mathrm{N}_{12}`$) is isoelectronic with $`\mathrm{C}_{60}^{+12}`$ ($`\mathrm{C}_{60}^{12}`$). Replacing 12 carbon atoms by 12 boron atoms results in removing 10 electrons from the HOMO ($`h_u`$ symmetry, $`E_h=5.99`$ eV) of $`\mathrm{C}_{60}`$ and 2 electrons from the 2nd HOMO ($`g_g`$ symmetry) of $`\mathrm{C}_{60}`$. In contrast, replacing 12 carbon atoms by 12 nitrogen atoms leads to a complete filling of 6 electrons in the LUMO ($`t_{1u}`$ symmetry, $`E_l=3.22`$ eV) of $`\mathrm{C}_{60}`$ and 6 electrons in the 2nd LUMO ($`t_{1g}`$ symmetry) of $`\mathrm{C}_{60}`$. Quantitatively, Fig.2 shows the orbital energies of $`\mathrm{C}_{48}\mathrm{B}_{12}`$ and $`\mathrm{C}_{48}\mathrm{N}_{12}`$. They are different from each other because of the difference in valency of B and N atoms. For $`\mathrm{C}_{48}\mathrm{N}_{12}`$, the HOMO is a doubly-degenerate level of $`\mathrm{a}_\mathrm{g}`$ symmetry, while the LUMO is a non-degenerate level with $`\mathrm{a}_\mathrm{u}`$ symmetry. For $`\mathrm{C}_{48}\mathrm{B}_{12}`$, the HOMO and LUMO are doubly-degenerate levels of $`\mathrm{a}_\mathrm{u}`$ and $`\mathrm{a}_\mathrm{g}`$ symmetries, respectively. These calculations show that $`\mathrm{C}_{48}\mathrm{B}_{12}`$’s LUMO is just above $`\mathrm{C}_{48}\mathrm{N}_{12}`$’s HOMO and the approximate potential difference is about 1.63 eV. Thus, both $`\mathrm{C}_{48}\mathrm{B}_{12}`$ and $`\mathrm{C}_{48}\mathrm{N}_{12}`$ have precisely the properties required for a proper rectifier behavior as will be demonstrated below.

Various experiments have shown that individual molecules can act as molecular rectifiers . The inset in Fig.3(a) shows a rectifier molecule consisting of a donor $`\mathrm{C}_{48}\mathrm{N}_{12}`$ and a acceptor $`\mathrm{C}_{48}\mathrm{B}_{12}`$ connected by a tunneling bridge (a $`\sigma `$-electron system, $`(\mathrm{CH}_2)_6`$). The bridge forms a potential barrier that isolates the donor from the acceptor on the time-scale of electron motion to or from the electrodes (see review  about synthesis techniques which have successfully bridged $`\mathrm{C}_{60}`$ dimers through electroactive spacers). Following the Aviram-Ratner scheme, we expect that electron current would pass from the cathode to the anode via acceptor, bridge and donor. Since $`\mathrm{C}_{48}\mathrm{B}_{12}`$’s LUMO is -4.2 eV, we may choose metallic carbon nanotubes (work function $`\varphi 5`$ eV), gold, silver or copper ( $`\varphi 4.4,4.7,4.8`$ eV, respectively ) as the electrodes. As discussed before, $`\mathrm{C}_{48}\mathrm{B}_{12}`$ and $`\mathrm{C}_{48}\mathrm{N}_{12}`$ have the properties required for proper rectifier behavior. An initial B3LYP/3-21G calculation for the full, connected fullerene/spacer/fullerene system shows that the covalent bond between the spacer and doped fullerenes shifts the electronic structure of the doped fullerenes. However, the HOMO and LUMO of the accepter-spacer-donor system are still localized on the acceptor and donor sides, respectively, ensuring the desired rectifier behavior. The same properties hold for the acceptor/bridge/donor system. Thus, a smaller threshold voltage for conduction is expected in one direction than in the other direction. As an example, we have investigated the rectifying characteristic of the $`\mathrm{C}_{48}\mathrm{B}_{12}/\mathrm{C}_6\mathrm{H}_{14}/\mathrm{C}_{48}\mathrm{N}_{12}`$ system by using non-equilibrium Green function theory in conjunction with a density functional based tight-binding model. The rectifier molecule is connected to gold contacts that we treat in the $`s`$-wave approximation with a constant density of states near to the Fermi level. The calculated current/applied potential characteristic is shown in Fig.3(a). A typical rectification characteristic is obtained with a turn-on bias close to the energy difference between $`\mathrm{C}_{48}\mathrm{B}_{12}`$’s LUMO and $`\mathrm{C}_{48}\mathrm{N}_{12}`$’s HOMO. Recently, Joachim et al. have successfully observed the electrical current flowing through an individual $`\mathrm{C}_{60}`$ molecule with a scanning tunneling microscope (STM). In light of this progress and the calculated size of $`\mathrm{C}_{48}\mathrm{X}_{12}`$ (X=N,B), our rectifier molecule also should be identifiable by STM.

Heterojunctions for molecular electronics can also be formed using $`\mathrm{C}_{60n}\mathrm{X}_n`$ donors and acceptors. Hybrid nanostructures formed by filling a SWNT with $`\mathrm{C}_{60}`$ have been observed, and can be superconducting  or metallic. Our DFT calculations, done within the local density approximation with double numerical basis including d-polarization function, show that the incorporation of the C48B12 or C48N12 into a (10,10)  or (17,0) SWNT is energetically favorable. About 2.4 eV binding energy per molecule is gained after the C48B12 or C48N12 molecule is inserted periodically inside the SWNTs. In our calculations of doped SWNTs with a one-dimensional periodic boundary condition along the tube axis, three unit cells of a (10,10) tube or two unit cells of a (17,0) tube are included in the supercell with one $`\mathrm{C}_{48}\mathrm{B}_{12}`$ or $`\mathrm{C}_{48}\mathrm{N}_{12}`$. Charge analysis found that placing an acceptor $`\mathrm{C}_{48}\mathrm{B}_{12}`$ into a (17,0) tube puts a +0.67 $`|e|`$ charge on the SWNT, while incorporating a donor $`\mathrm{C}_{48}\mathrm{N}_{12}`$ into a (17,0) SWNT puts a -0.39 $`|e|`$ charge on the SWNT. Similar results are obtained for a (10,10) SWNT. Hence, putting $`\mathrm{C}_{48}\mathrm{B}_{12}`$ into a semiconducting tube results in a $`p`$-type region on the SWNT, while filling of donors $`\mathrm{C}_{48}\mathrm{N}_{12}`$ into a semiconducting tube leads to a $`n`$-type region on the SWNT. Thus, it is possible to use them to get $`n`$-$`p`$-$`n`$ and $`p`$-$`n`$-$`p`$ transistors. As shown in Fig.3(b), doping a (17,0) semiconducting SWNT with $`\mathrm{C}_{48}\mathrm{B}_{12}`$ and $`\mathrm{C}_{48}\mathrm{N}_{12}`$ should also make a SWNT-based $`p`$-$`n`$ junction.

FIG.3: (a) Calculated current through a rectifier molecule (inset) which consists of $`\mathrm{C}_{48}\mathrm{B}_{12}`$ (left, red for B atoms) and $`\mathrm{C}_{48}\mathrm{N}_{12}`$ (right, blue for N atoms) connected by $`\mathrm{C}_6\mathrm{H}_{14}`$ (middle, green for H atoms) and is coupled to two Au electrodes via S atoms on the two sides. The two C atoms in the ends of the $`\mathrm{C}_6\mathrm{H}_{14}`$ bridge are connected to the N (B) atoms in $`\mathrm{C}_{48}\mathrm{N}_{12}`$ ($`\mathrm{C}_{48}\mathrm{B}_{12}`$) respectively. A B3LYP/3-21G calculation shows that the connected rectifier molecular superstructure is stable relative to the separated, individual components by about 19.545 kcal/mol energy. (b) a prototype for $`\mathrm{C}_{48}\mathrm{X}_{12}`$@(17,0) SWNT-based (X=B,N) $`p`$-$`n`$ junction.

In conclusion, we show that acceptor $`\mathrm{C}_{60n}\mathrm{B}_n`$ and donor $`\mathrm{C}_{60m}\mathrm{N}_m`$ pairs needed for molecular electronics can be obtained by properly controlling the number $`n`$ and $`m`$ of the substitutional dopants in $`\mathrm{C}_{60}`$. We demonstrated the rectifying characteristic of a molecular rectifier built from our engineered acceptor/donor pairs. Heterojunctions for molecular electronics can be made by inserting these dopants into semiconducting carbon nanotubes. Efficiently synthesizing those acceptor/donor pairs would be of great experimental interest within reach of today’s technology. Very recently, small nitrogen-substitutionally-doped fullerenes have been reported,  showing experimental progress in this direction.

One of us (R.H.X.) thanks the HPCVL at Queen’s University for the use of its parallel supercomputing facilities. VHS gratefully acknowledges support from NSERC. ADC acknowledges support from the European DIODE Network and the Italian MURST-CNR under the project “Nanoelectronics”.