text:
# Modeling Nonlinear Dynamical Systems with Delay-differential Equations.

## 1 Introduction

How the information about the outside world is stored in the brain is still in many ways an open question. Clearly, representation of the totality of the information involves some sort of modeling of the environment by the brain under the specific conditions of its operation, e.g., the limited number of the neurons and synapses, the inherent delays in signal propagation along the neural fibers, etc. Assuming that the environment is described by some nonlinear deterministic dynamical system and that the brain is described by another dynamical system, this question can be rephrased as the question of how one dynamical system can model another. Since the environment is typically much more complicated than the brain, the modelling task should be impossible simply on the count of the difference in the number of the degrees of freedom for the two dynamical systems. A possible solution to this puzzle is to assume that the environment is actually a collection of weakly interacting subsystems, each of which has a lower complexity than the brain. A good strategy to learn complicated environment would then be to model the strongly coupled degrees of freedom group by group. Decomposition of the environment into a collection of components that interact weakly presumably should involve some version of nonlinear factor analysis extended to deterministic systems.

Even if the complexity of the environment is less than that of the system that models it a fundamental question arises about in what sense and how one dynamical system can model another. In addition one can ask whether there exist a class of dynamical systems that in some sense are universal in their modeling properties in the sense that a large class of nonlinear dynamical systems that represent possible environments can be modeled using essentially the same system. Of course, in order for this to make sense in explaining the properties of the brain, the modeling has to be physically realizable.

Leaving investigation of the deterministic nonlinear factor analysis to another publication , in this paper we present arguments that when the dimension of the environment is less or equal to that of the modeling system such universal modeling systems exist and that it is by the essential use of the delays that the universality can be realized. We prove that under certain conditions any nonlinear dynamical system can be modeled by a linear delay-differential dynamical system with constant coefficients, in the sense that any finite set of trajectories of the environment can be loaded into periodic solutions of the delay-differential equation. The loading turns out to be equivalent to solving an essentially linear problem and is achieved by adjusting the coefficients of the delay-differential equation in a way that is plausible from the biological point of view. With the biological applications in mind we also extend some of the results to a class of nonlinear systems with delays.

Modeling the environment is usually referred to as construction memory models of the environment. Since the early eighties much work has been done on memory models under the assumption that time-independent memories can be represented as fixed points of systems of nonlinear differential equations (ODE) that possess a Liapunov function and with the evolution given by

$$\dot{x}\left(t\right)=F\left(x\right),x(t)R^N$$
(1)
for some function $`F\left(x\right)`$. Once a Liapunov function is given, the fixed points of the evolution can be identified with the $`N`$-dimensional stored patterns. Domains of attraction of the fixed points then can be viewed as all the patterns that will be ”associatively” recalled dynamically as the stored pattern represented by the fixed points. Especially useful within this approach proved to be the Hopfield model  described as a dynamical system with evolution

$$\dot{x}\left(t\right)=E_\mu x\left(t\right)+A\sigma (x\left(t\right))+y,$$
(2)
where $`E_\mu =diag\{\mu _i\},\mu _i>0`$ is an $`N\times N`$ diagonal matrix, $`A`$ \- is an $`N\times N`$ matrix, called the matrix of weights, $`\sigma :x_i\sigma _i(x_i)`$ is a component-wise nonlinear transformation and the constant vector $`yR^N`$ describes time-independent external inputs to the system. If $`\sigma ^{}(z)0`$ and matrix $`A`$ is symmetric, or if the symmetric part of $`A`$ is non-positive definite , then the system $`\left(\text{2}\right)`$ has a Liapunov function $`V\left(x\right)`$ defined on the configuration space such that on the trajectories $`\dot{V}(x(t))<0,`$ leading to the convergence of all trajectories to a set of fixed points. An extensive theory exists about the properties of memory storage, especially for random patterns .

A certain amount work has been done for the systems with delays. Existence of a Liapunov function can be proven for the Hopfield model with delays if the delays are not too large . However, the fixed-points paradigm is difficult to apply to the time dependent patterns. Time-variable patterns, which in our approach are represented by the trajectories of nonlinear dynamical systems, are ubiquitous in the environment. Therefore, understanding the underlying principles of their storage is important when constructing plausible models for perception and memory. Some interesting effects of the presence of delays on the structure of the attractors leading to multistability in networks of two Hodgkin-Huxley type spiking neurons were discussed heuristically and numerically in . From our point of view the model that is considered there is that of a system of two linear delay differential equations with time variable coefficients. In this paper we do not restrict ourselves to the models of spiking neurons but consider arbitrarily large systems of linear and nonlinear delay-differential equations with slowly varying coefficients. Systems with time variable coefficients will be considered in detail elsewhere. Nevertheless, our results are in general agreement with the numerical simulations in  and provide a new perspective on their interpretation.

Here we present an approach to the time-dependent pattern storage that describes sets of interacting ”neurons” as systems of delay-differential equations (DDE) with multiple delays. Since delays in nerve impulse propagations are common, the approach is more plausible biologically than, for example, the Hopfield model where instantaneous communication among the neurons is assumed. Mathematically the delays appear after integrating out some internal degrees of freedom in systems of Hodgkin-Huxley equations . Although more complicated to analyze, the DDE have some properties that ordinary differential equations do not have. For example, for a solution of a DDE to be unique, one has to specify an initial condition which is a function defined on a finite time range . This is a key feature of DDE that we exploit.

Our starting point is to include explicitly the multiple delays $`\tau _1,\mathrm{},\tau _L`$ in the dynamics of memory model so that instead of $`\left(\text{1}\right)`$ we obtain a system of nonlinear DDE

$$\dot{x}\left(t\right)=F(x\left(t\right),x\left(t\tau _1\right),\mathrm{},x\left(t\tau _L\right)).$$
(3)
One example of such a system is a generalization of the Hopfield model to include delays. This was considered in

$$\dot{x}\left(t\right)=\underset{k=0}{\overset{L}{}}A_k\sigma \left(x\left(t\tau _k\right)\right)+y(t),$$
(4)
where as before the $`A_k`$ and are $`N\times N`$ weight matrices with $`A_0=E_\mu ,`$ and the delays are defined so that $`\tau _0=0,\tau _k`$ $``$ $`\tau _{k+1}`$. In contrast to $`\left(\text{2}\right)`$ we assume the external inputs $`y(t)`$ as time-variable. In models describing networks of oscillating neurons one can also consider $`\tau _k`$ as an $`N\times N`$ matrix with each element $`\left(\tau _k\right)_{ij}`$ describing the delay in signal propagation from neuron $`\left(j\right)`$ to neuron $`\left(i\right)`$ via $`\mathrm{"}k\mathrm{"}`$th pathway. However, after appropriate redefinition of matrices $`A_k`$ and nonlinearities $`\sigma \left(z\right)`$ one can recover $`\left(\text{4}\right)`$ at the expense of increasing $`L`$ to the number of different values of all matrix elements $`\tau _k`$ of all $`\left(L+1\right)`$ matrices. We assume that $`y(t)`$ is a trajectory of some nonlinear dynamical system with local evolution given by $`\dot{y}(t)=G(y).`$ The detailed nature of $`G(y)`$ is not important for the discussion. When $`\sigma _i(z)=z`$ we obtain a system of linear DDE. Writing it in components we obtain

$$\dot{x}_i=\underset{k=0}{\overset{L}{}}\left(A_k\right)_{ij}x_j\left(t\tau _k\right)+y_i(t),i=1,\mathrm{},N.$$
(5)
The main point of this paper is to prove the existence of the periodic solutions of systems of homogeneous linear DDE that are periodic extensions of their initial conditions. That is, we show that, given a set of fixed delays $`\tau _k`$, the weight matrices $`A_k`$ can be chosen in such a way that the initial condition when extended periodically to $`t>\tau _L`$ is the solution of $`\left(\text{5}\right)`$. With appropriate choice of parameters the choice of $`A_k`$ is unique. We also show that if the inhomogeneous part of a DDE is given by a trajectory of a nonlinear dynamical system, the trajectory can be loaded approximately, but with arbitrarily small error into the initial condition for the DDE and hence stored as its periodic solution. With some modifications these results can be carried over to a class of nonlinear DDE. Using the method we describe, with a sufficiently large $`N,`$ any finite set of trajectories of a nonlinear dynamical system, for example any finite set of periodic orbits, can be encoded as periodic solutions of a linear DDE. Since many dynamical systems, are uniquely characterized by the sets of their periodic orbits, we can speak of linear modelling of nonlinear systems.

The paper is structured as follows. In the next section we give a brief overview of the differences between the linear ODE and the linear DDE from mathematical point of view. Section 3 contains discussion of loading the trajectories into the initial conditions. In Section 4 we construct solutions for the periodic extension problem. In Section 5 we discuss the adaptation dynamics for the weight matrices $`A_k`$. Section 6 is a summary.

## 2 Delay-differential Equations

To specify a solution of a differential equation uniquely one needs to pose the initial value problem. There is a big difference in how the initial value problem is posed for ODE and for DDE. For the ODE the initial value problem is given by

$`x(\varphi ,y)\left(t\right)`$ is a solution of $`\left(\text{1}\right)`$ for $`t>0;`$

$`x(\varphi ,y)\left(0\right)=\varphi R^N.`$

For the DDE of retarded type $`\left(\text{3}\right)`$ with the maximum delay $`\tau _L`$ the initial value problem is formulated as

$`x(\varphi ,y)\left(t\right)`$ is a solution of $`\left(\text{3}\right)`$ for $`t>\tau _L;`$

$`x(\varphi ,y)\left(0\right)=\varphi \left(t\right)`$ $`R^N`$ for $`t[0,\tau _L].`$

Hence to specify a unique solution of a DDE one needs an initial value from a space of functions defined on the interval $`[0,\tau _L]`$. Exploiting the difference is one key feature of our approach.

Let us consider in more detail the difference in the linear case. Similar to homogeneous linear ODE we can search for the exponential solutions of homogeneous DDE. As for ODE these correspond to zeroes of the characteristic equation that is obtained using the Laplace transform. For linear ODE the characteristic equation is polynomial and has exactly $`N`$ complex roots. The situation is different for linear DDE. The characteristic equation becomes transcendental and is given by

$`det\left(H\left(\lambda \right)\right)`$ $`=0,`$ (6)
$`H\left(\lambda \right)`$ $`=\lambda I{\displaystyle \underset{k=0}{\overset{L}{}}}A_k\mathrm{exp}\left(\lambda \tau _k\right).`$
There are generically infinite number of complex roots $`\lambda _p`$ of $`\left(\text{6}\right).`$ Let $`N(a,b),`$ $`a,bR,`$ $`ab<\mathrm{},`$ be the number of the complex roots of $`\left(\text{6}\right)`$ in the vertical strip of the complex plain with real part of each point contained in $`[a,b].`$ Then, some of useful properties of the roots are

(i) Re$`\left(\lambda _p\right)c<\mathrm{}`$ if all $`\tau _k0,`$ for all $`k,i,j;`$

(ii) Re$`\left(\lambda _p\right)\mathrm{}`$ $`p\mathrm{};`$

(iii) $`N(a,b)<\mathrm{}`$.

Solutions of linear inhomogeneous DDE can be conveniently written down in an integral form using the notion of the fundamental solution. The fundamental solution is defined as the unique matrix solution of $`\left(\text{5}\right)`$ with $`y\left(t\right)=0`$ and with the following initial condition

$`X\left(t\right)=0`$ $`t[0,\tau _L];`$

$`X\left(t\right)=I`$ $`t=\tau _L,`$ where $`I`$ is the unit matrix.

Using inverse Laplace transform, the fundamental solution can be expressed in terms of the characteristic (matrix) function $`H\left(\lambda \right)`$

$$X\left(t\right)=\frac{1}{2\pi i}_{\left(c\right)}𝑑\lambda e^{\lambda t}H^1\left(\lambda \right),$$
(7)
where $`\left(c\right)`$ is the contour $`(ci\mathrm{},c+i\mathrm{})`$ in the complex plane such that $`\mathrm{Re}\left(\lambda _p\right)<c`$ and where the Laplace transform and its inverse are defined by

$`\widehat{x}\left(\lambda \right)`$ $`={\displaystyle _0^{\mathrm{}}}𝑑\lambda \mathrm{exp}\left(\lambda t\right)x\left(t\right),`$ (8)
$`x\left(t\right)`$ $`={\displaystyle \frac{1}{2\pi i}}{\displaystyle _{\left(c\right)}}𝑑\lambda \mathrm{exp}\left(\lambda t\right)\widehat{x}\left(\lambda \right).`$ (9)
With these definitions the unique solution of the initial value problem for $`\left(\text{5}\right)`$ can be written as

$`x(\varphi ,y)\left(t\right)`$ $`=X\left(t\tau _L\right)\varphi \left(\tau _L\right)+{\displaystyle _{\tau _L}^t}𝑑\tau X\left(t\tau \right)y\left(\tau \right)`$ (10)
$`+{\displaystyle \underset{k=1}{\overset{L}{}}}{\displaystyle _{\tau _L\tau _k}^{\tau _L}}𝑑\tau X\left(t\tau \tau _k\right)A_k\varphi \left(\tau \right).`$
If we put $`L=0`$ then the third term in $`\left(\text{10}\right)`$ vanishes and we recover the solution of a linear inhomogeneous ODE in terms of its fundamental solution. Note that in $`\left(\text{10}\right)`$ the inhomogeneous part $`y\left(t\right)`$ and the initial condition $`\varphi \left(t\right)`$ play a similar role. We shall exploit the similarity below when we generate initial conditions from the inhomogeneous part of the equation.

## 3 Loading Trajectories into the Initial Conditions

In this section we describe how to generate the initial condition for a DDE from its inhomogeneous part. The main idea is to consider the inhomogeneous part to be nonzero only on $`[0,\tau _L]`$ and generate the initial condition recursively by adding more and more terms to the equation successfully.

Let us first consider the linear DDE. We divide the interval $`[0,\tau _L]`$ into adjoining intervals $`\left[\tau _{k,}\tau _{k+1}\right],k=0,\mathrm{},L1`$ and for time $`t<`$ $`\tau _m`$ truncate the full DDE to

$$\dot{x}\left(t\right)=\underset{k=0}{\overset{m}{}}A_kx\left(t\tau _k\right).$$
(11)
This can be done under the assumption that the signals with larger delays did not arrive yet at the ”neuron” to influence its state.

When $`m=0`$ $`\left(\text{11}\right)`$ reduces to a linear ODE which can be integrated for $`t(0,\tau _1]`$ using the initial condition $`\varphi _0`$ at $`t=0`$ and the values of the inhomogeneous part $`y(t)`$. The result of the integration can be considered as the initial condition the truncated DDE defined on $`t(\tau _1,\tau _2]`$, which enables us to integrate the DDE with one delay $`\tau _1`$on the segment $`\left[\tau _{1,}\tau _2\right]`$ using $`y(t)`$ defined on $`\left[\tau _{1,}\tau _2\right]`$ only. Proceeding in this fashion step by step we can construct the initial condition on the entire segment $`[0,\tau _L]`$. Taking the inhomogeneous part $`y(t)`$ to be zero outside of $`[0,\tau _L]`$ we can interpret the result of this procedure as a linear mapping of trajectories into the initial conditions, $`I_L:`$ $`y(t)`$ $`\mathrm{\Phi }\left(t\right),t`$ $`[0,\tau _L]`$.

To describe the iteration procedure in more detail we give the initial iteration step followed by the $`\mathrm{"}m\mathrm{"}`$th iteration. For the initial step we need to generate the initial condition for

$`\dot{x}\left(t\right)`$ $`=A_0x\left(t\right)+A_1x\left(t\tau _1\right)+y(t),`$ (12)
$`y(t)`$ $`=0,t[0,\tau _1].`$
These are given by the solution of the initial value problem

$`\dot{x}\left(t\right)`$ $`=A_0x\left(t\right)+y(t),`$ (13)
$`x(0)`$ $`=\varphi _0.`$
Using the fundamental solution for this equation the solution for $`t>0`$ and the initial condition for $`\left(\text{12}\right)`$ can be written as

$`\mathrm{\Phi }_0\left(t\right)`$ $`=X_0\left(t\right)\varphi _0+{\displaystyle _0^t}𝑑sX_0\left(ts\right)y\left(s\right),`$ (14)
$`X_0\left(t\right)`$ $`=\mathrm{exp}\left(A_0t\right)`$
With this initial condition the solution for $`\left(\text{12}\right)`$ can be written in terms of its fundamental solution as

$`\mathrm{\Phi }_1\left(t\right)`$ $`=X_1\left(t\tau _1\right)\mathrm{\Phi }_0\left(\tau _1\right){\displaystyle _0^{\tau _1}}𝑑sX_1\left(ts\tau _1\right)A_1\mathrm{\Phi }_0\left(s\right)`$ (15)
$`+{\displaystyle _{\tau _1}^t}𝑑sX_1\left(ts\right)y\left(s\right).`$
$`\mathrm{\Phi }_1\left(t\right)`$ can be considered as the initial condition for the 2nd step of iteration, the result of which is the initial condition defined on the segment $`\left[\tau _{1,}\tau _2\right].`$ The general form of iteration is easily deduced. Namely,

$`\mathrm{\Phi }_k\left(t\right)`$ $`=X_k\left(t\tau _k\right)\mathrm{\Phi }_{k1}\left(\tau _k\right){\displaystyle _{\tau _{k1}}^{\tau _k}}𝑑sX_k\left(ts\tau _k\right)A_k\mathrm{\Phi }_{k1}\left(s\right)`$ (16)
$`+{\displaystyle _{\tau _k}^t}𝑑sX_k\left(ts\right)y\left(s\right).`$
with the final step resulting in the initial condition

$$\mathrm{\Phi }_L\left(t\right)=\{\mathrm{\Phi }_k\left(t\right),t[\tau _k,\tau _{k+1}],k=0,\mathrm{},L1\}.$$
(17)
Next note that to enable its periodic extension the initial condition has to be periodic, since it has to be defined on a closed segment. Hence we obtain a constraint of the form

$$\mathrm{\Phi }_L\left(\tau _L\right)=\mathrm{\Phi }_L\left(0\right)\varphi _0.$$
(18)
This condition can be used to eliminate an unknown parameter $`\varphi _0`$ and make trajectory loading unique. For example, when there is only one delay we obtain

$$\varphi _0=\left(\mathrm{exp}\left(A_0\tau _L\right)1\right)^1_0^{\tau _L}𝑑s\mathrm{exp}\left(A_0s\right)y(s).$$
(19)
When multiple delays are present $`\varphi _0`$ can be determined similarly. The null space of the mapping $`I_L`$ is also of interest. It describes all the trajectories that cannot be loaded into the periodic extensions of the initial conditions. For the example with only one delay this space consists of all functions with the Laplace transform satisfying

$$\widehat{y}\left(\lambda \right)+\varphi _0=0$$
(20)
For arbitrary set of delays the null space still forms a one-parameter family and, therefore, its existence excludes only a very small set of all trajectories from loading.

Obviously, the detailed loading procedure described above works only for the linear DDE. However the same principle can be applied to the nonlinear DDE as well. If the delays are known, then one can proceed with the same iterative scheme by numerical integration. We speculate that some sort of the analog integration version of the procedure above might be used in the brain to implement the loading. With this in mind we can interpret the time $`\tau _L`$ as the ”attention span” of the population of the neurons that encode a particular time-dependent pattern. For neurons with a large number of synapses and, hence a large number of random delays, a continuous version of the iteration in $`\left(\text{16}\right)`$ can be easily written.

## 4 Periodic Extensions of the Initial Conditions

Having shown how trajectories can be loaded into the initial conditions of the linear DDE, we now construct solutions of the DDE that are the periodic extension of their initial conditions. First we consider the linear case and then discuss the modifications added by the presence of nonlinearities.

### 4.1 Linear DDEs

Consider a Fourier series representation of a space of periodic and continuous initial conditions on $`[0,\tau _L]`$

$`\varphi \left(t\right)`$ $`={\displaystyle \underset{n=Q}{\overset{Q}{}}}\mathrm{exp}\left(i\rho _nt\right)b_n,t[0,\tau _L],`$ (21)
$`\rho _n`$ $`={\displaystyle \frac{2\pi }{\tau _L}}n,`$
where $`b_n`$ are fixed $`N`$-dimensional complex-valued vectors such that under complex conjugation$`\overline{b}_n=`$ $`b_n`$ (to ensure that $`\varphi \left(t\right)`$ is real-valued). If the initial condition does not belong to this space we shall consider $`\left(\text{21}\right)`$ as an approximation of the initial condition with the error of approximation determined by the truncation of the infinite Fourier series that represents it at the ”$`Q\mathrm{"}`$th term of the expansion.

Let us extend the domain of definition of the initial condition $`\varphi \left(t\right)`$ to all $`t0`$ by treating $`\varphi \left(t\right)`$ as a continuous periodic function for $`t0`$ with $`\varphi \left(t+\tau _L\right)=\varphi \left(t\right).`$ Assume now that all delays $`\tau _k`$ are fixed but we are free to vary the matrices $`A_k.`$ The problem of finding periodic extensions is then to find such $`A_k`$ that for $`t0`$ the function $`\varphi \left(t\right)`$ is a solution of the homogeneous equation

$$\dot{\varphi }\left(t\right)=\underset{k=0}{\overset{L}{}}A_k\varphi \left(t\tau _k\right).$$
(22)
After the substitution of $`\left(\text{21}\right)`$ into $`\left(\text{22}\right)`$ we obtain a set of linear equations on matrices $`A_k`$

$$\left(i\rho _nI\underset{k=0}{\overset{L}{}}A_k\mathrm{exp}\left(i\rho _n\tau _k\right)\right)b_n=0,n=0,\mathrm{},Q.$$
(23)
Note that the equations for $`n=Q,\mathrm{},1`$ are redundant, since $`A_k`$ are real and hence the additional equations can be obtained from $`\left(\text{23}\right)`$ by complex conjugation. Altogether, because the $`n=0`$ equation is real, $`\left(\text{23}\right)`$ is equivalent to $`\left(2Q+1\right)N`$ real equations on $`\left(L+1\right)N^2`$ elements of $`A_k,`$ provided that all elements of $`A_k`$ are generically non-zero. In the case a delay $`\tau _k`$ is given by a random value of delays in propagation from $`\mathrm{"}i\mathrm{"}`$th to $`\mathrm{"}j\mathrm{"}`$ neurons, with probability one only one of the elements of $`A_k`$ is non-zero (see the discussion in the Introduction). However, in such a case, if $`L^{}+1`$ is the number of $`N\times N`$ delay matrices, then $`L+1=\left(L^{}+1\right)N^2`$ and, therefore, the count of unknowns is the same. We conclude that $`\left(\text{23}\right)`$ can have a unique solution for $`A_{k\text{ }}`$only when $`\left(2Q+1\right)=\left(L+1\right)N.`$ Since $`Q`$ determines the number of terms in the Fourier expansion and, therefore, the accuracy of representation of an arbitrary initial condition, this relation means that for a given number $`L`$ of delays and $`N`$ the number of neurons involved in modeling one can achieve only limited accuracy of representation. Namely, the number of expansion terms $`Q`$ is bounded so that $`\left(2Q+1\right)\left(L+1\right)N.`$

The system $`\left(\text{23}\right)`$ can be solved by noticing that it implies that vectors $`b_n`$ are eigenvectors of linear combinations of $`A_k`$ with eigenvalues $`i\rho _n.`$ Unlike the typical eigenvalue problem where $`A_k`$ are known and one needs to find $`\rho _n`$ and $`b_n`$, here the situation is reversed: one needs to find $`A_k`$ assuming that $`\rho _n`$ and $`b_n`$ are known. If an eigenvalue $`\lambda `$ and an eigenvector $`b`$ are given then all the matrices $`D`$ for which they solve the eigenvalue problem can be written as

$$D=\lambda I+B\left(Ib^2b\overline{b}\right),$$

where $`B`$ is arbitrary $`N\times N`$ matrix and $`b^2=_{i=1}^N\left|b_i\right|^2`$ and $`b\overline{b}`$ is the matrix that is the outer product of $`b`$ with its complex conjugate. For given eigenvalue and eigenvector matrix $`D`$ has $`N^2N`$ free parameters. Therefore, an equivalent way to write $`\left(\text{23}\right)`$ is

$`{\displaystyle \underset{k=0}{\overset{L}{}}}R_{nk}A_k`$ $`=i\rho _nI+B_n\left(Ib_n^2b_n\overline{b}_n\right),n=0,\mathrm{},Q,`$ (24)
$`R_{nk}`$ $`=\mathrm{exp}\left(i\rho _n\tau _k\right),`$ (25)
where $`B_n`$ is an $`N\times N`$ complex matrix, which effectively has $`N^2N`$ free parameters. This linear system of equations always has solutions for appropriate choice of parameters $`Q,L,N.`$ For some choices the solution is unique. This concludes the proof of existence of the periodic extensions.

Let us consider $`(\text{23},\text{24})`$ for various choices of parameters. When $`N=1`$ the second term in the RHS of $`\left(\text{24}\right)`$ vanishes and the dependence on $`b_n`$ drops out. Consequently, the system $`\left(\text{23}\right)`$ becomes a system of $`\left(2Q+1\right)`$ real linear equations on $`L+1`$ real parameters $`A_k`$

$`{\displaystyle \underset{k=0}{\overset{L}{}}}S_{nk}A_k`$ $`=\rho _n,n=1,\mathrm{},Q,`$ (26)
$`{\displaystyle \underset{k=0}{\overset{L}{}}}C_{nk}A_k`$ $`=0,n=0,1,\mathrm{},Q,`$ (27)
$`S_{nk}`$ $`=\mathrm{sin}\left(\rho _n\tau _k\right),`$ (28)
$`C_{nk}`$ $`=\mathrm{cos}\left(\rho _n\tau _k\right).`$ (29)
When $`L=2Q,`$ provided $`detT0,`$ the equations can be solved uniquely by inversion of the matrix $`T`$

$$T=\left(\begin{array}{ccc}\mathrm{sin}\left(\rho _1\tau _0\right)& \mathrm{}& \mathrm{sin}\left(\rho _1\tau _L\right)\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{sin}\left(\rho _Q\tau _0\right)& \mathrm{}& \mathrm{sin}\left(\rho _Q\tau _L\right)\\ 1& \mathrm{}& 1\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{cos}\left(\rho _Q\tau _0\right)& \mathrm{}& \mathrm{cos}\left(\rho _Q\tau _L\right)\end{array}\right)$$
(30)
For arbitrary $`N`$ , $`L=2Q,`$ and $`detT0`$ we can write down the solutions for $`\left(\text{23}\right)`$ $`R`$ as

$`A_k`$ $`={\displaystyle \underset{k=0}{\overset{2Q}{}}}T_{kn}^1\mathrm{\Gamma }_n,k=0,1,\mathrm{},L.`$ (31)
$`\mathrm{\Gamma }_n`$ $`=\left(\rho _nI+\mathrm{Im}\left(B_n\left(Ib_n^2b_n\overline{b}_n\right)\right)\right),n=1,\mathrm{},Q`$ (32)
$`\mathrm{\Gamma }_n`$ $`=\mathrm{Re}\left(B_n\left(Ib_n^2b_n\overline{b}_n\right)\right),n=0,1,\mathrm{},Q`$ (33)
Since the matrices $`\mathrm{\Gamma }_n`$ contain $`\left(N^2N\right)`$ free parameters, to obtain a unique solution one needs to increase the value of $`Q`$ by factor of $`\left(N^2N\right)`$ to constrain the additional degrees of freedom.

Some remarks about the meaning of the periodic extensions should be made. Since the initial conditions were taken as a linear combination of exponential solutions, the existence of the extensions is equivalent to the existence of $`2Q+1`$ pair-wise conjugate zeroes of the characteristic function that are located on the imaginary axis. As mentioned in the previous section, there can be only a finite number of zeroes of the characteristic function in any finite width vertical strip of the complex plane. Hence, there could be only a finite number of zeros of the characteristic function lying on the imaginary axis. As a result, $`Qmust`$ be finite and only finite-dimensional spaces of initial conditions can be extended periodically. A practical consequence of this is that one can store a trajectory as a periodic extension only approximately, with the error of approximation given by the error induced by truncating the Fourier series. However, as follows from the theory of Fourier expansions, the error can be made arbitrarily small at least in $`L^2`$ norm by increasing $`N.`$

#### 4.1.1 Non-linear DDE

Let us now consider the existence of the periodic extensions in the presence of nonlinearities. Take, for example, the extended Hopfield model with dynamics

$$\dot{x}=\underset{k=0}{\overset{L}{}}A_k\sigma \left(x\left(t\tau _k\right)\right)+y(t)$$
(34)
Proceeding as before we expand the periodic the initial condition in the Fourier series for $`t[0,\tau _L]`$ as in $`\left(\text{21}\right)`$ and obtain that for $`\varphi \left(t\right)`$ to be a solution of homogeneous DDE it needs to satisfy

$$\dot{\varphi }\left(t\right)=\underset{k=0}{\overset{L}{}}A_k\sigma \left(\varphi \left(t\tau _k\right)\right)$$
(35)
Substitution of $`\left(\text{21}\right)`$ into $`\left(\text{35}\right)`$ and expansion of the nonlinear term in Fourier series yields a system of equations on the unknown matrices $`A_k`$

$$\left(i\rho _nb_n\underset{k=0}{\overset{L}{}}A_k\mathrm{exp}\left(i\rho _n\tau _k\right)\mathrm{\Lambda }_n(\tau _k;\left\{b_m\right\})\right)=0,n=Q_1,\mathrm{},Q_1$$
(36)
where $`Q_1`$ is not necessarily equal to $`Q`$ and for each $`p=1,\mathrm{}.N`$ the coefficient $`\left(\mathrm{\Lambda }_n(\tau _k;\left\{b_m\right\})\right)_p`$ is defined by

$$\mathrm{\Lambda }_n(\tau _k;\left\{b_m\right\})=\frac{1}{\tau _L}_0^{\tau _L}𝑑t\mathrm{exp}\left(i\rho _n\left(t\tau _k\right)\right)\sigma \left(\underset{m=Q}{\overset{Q}{}}\mathrm{exp}\left(i\rho _m\left(t\tau _k\right)\right)b_m\right)$$
(37)
The additional time-independent factor $`\mathrm{exp}\left(i\rho _n\tau _k\right)`$ in the definition of $`\mathrm{\Lambda }_n(\tau _k;\left\{b_m\right\})`$ ensures that if $`\sigma \left(z\right)=z`$ then $`\mathrm{\Lambda }_n(\tau _k;\left\{b_m\right\})=b_n.`$

Now the vector coefficients of the Fourier expansion $`\left(\text{21}\right)`$ cannot be factored out and the solution for $`A_k`$ will depend not only on the delays $`\tau _k`$ but also on $`b_n`$even for $`N=1`$. Such dependence is desirable in a memory model, since it implies that the choice of the parameters $`A_k`$ is pattern specific.

When $`Q_1`$ $``$ $`Q`$ and $`Q_1<\mathrm{}`$ we obtain essentially the same system of equations as before. One class of nonlinearities when this occurs is a set that may be called polynomial nonlinearities with

$$\sigma \left(z\right)=\underset{N_1}{\overset{N_2}{}}c_nz^n,0<N_1N_2<\mathrm{}$$
(38)
For example, if we assume that all trajectories that we wish to store are bounded, in order to introduce a sigmoid nonlinearity that is frequently used in memory models, one can take a cubic nonlinearity with

$$\sigma \left(z\right)=\left(1/3\right)\alpha ^2z^3+\beta ^2z,\alpha ,\beta >0$$
(39)
and load only such trajectories for which $`\sigma \left(z\right)\left(2/3\right)\left(\beta ^3/\alpha \right).`$ Substitution of $`\left(\text{39}\right)`$ into $`\left(\text{35}\right)`$ results in

$`\left(i\rho _nb_n{\displaystyle \underset{k=0}{\overset{L}{}}}A_k\mathrm{exp}\left(i\rho _n\tau _k\right)\left(\beta ^2b_n+{\displaystyle \frac{1}{3}}{\displaystyle \underset{p,q=Q}{\overset{Q}{}}}b_pb_qb_{npq}\right)\right)`$ $`=0,`$ (40)
$`n`$ $`=Q_1,\mathrm{},Q_1`$
As in the linear case this equation can be still considered as an eigenvalue problem for $`b_n`$, except for a nonlinear operator $`\stackrel{~}{A}=g\left(A\right).`$ On the other hand this equation can also be considered as a linear equation on $`A_k`$ for a given set of $`b_n.`$ One novel feature that appears in the nonlinear case is that since we have a cubic nonlinearity, in general, up to three different values for each $`b_n`$ are possible so that the term $`\left(\beta ^2b_n+\frac{1}{3}_{p,q=Q}^Qb_pb_qb_{npq}\right)`$ has the same value and hence the solution for $`A_k`$ is the same. This observation can be used for multiple memory storage of a set of $`M`$ desired memories $`\left\{b_n^B\right\},B=1,\mathrm{},M.`$ Consider the following system of equations

$$\left(\beta ^2b_n+\frac{1}{3}\underset{p,q=Q}{\overset{Q}{}}b_pb_qb_{npq}\right)=C_n,n=Q_1,\mathrm{},Q_1$$

This cubic system of $`2Q+1`$ equations has $`3^{2Q+1}`$ choices of solutions. As a result, in general, one can load $`3^{2Q+1}`$ memories corresponding to the same choice of matrices $`A_k.`$ The same argument applies to any polynomial nonlinearity defined by $`\left(\text{38}\right)`$ with the corresponding maximum number of memories growing as $`N_2^{2Q+1}.`$ Of course for this storage method to be practical a prescription should be given for a way to store a given set of memories for a given nonlinearity, rather then choosing a nonlinearity to fit the needed set of memories as we described above. In addition a thorough investigation of the basins of attraction of the multiple memories needs to be carried out. These and other questions shall be addressed in a future publication.

A problem can arise if $`Q_1`$ is sufficiently large, that is the system can become overdetermined if $`Q_1`$ is large enough. If $`Q_1`$ is finite the problem can be cured by choosing appropriately large $`L.`$ The real problem appears when nonlinearity is such that $`Q_1=\mathrm{}.`$ Then, one would expect that the system $`\left(\text{35}\right)`$ generically has no solutions. When $`Q_1=\mathrm{}`$ our approach can be nevertheless applied approximately in the following sense. If $`\mathrm{\Lambda }_n(\tau _k;\left\{b_m\right\})0`$ sufficiently fast as $`n\mathrm{},`$ we can truncate the system at a finite $`Q_{\mathrm{max}}`$ and the proceed as before with computing the $`A_k`$. Of course since we zeroed the coefficients $`\mathrm{\Lambda }_n(\tau _k;\left\{b_m\right\})`$ for $`n>Q_{\mathrm{max}}`$ in effect we substituted the original nonlinearity $`\sigma \left(z\right)`$ with another one $`\stackrel{~}{\sigma }\left(z\right)`$ and the periodic extensions obtained for $`\stackrel{~}{\sigma }\left(z\right)`$ will not be periodic solutions for the original nonlinearity. However if the original nonlinear DDE is stable with regard to small perturbations of the solutions, the aperiodicity will not grow in time and will remain small. Therefore we still be able to store time-variable patterns, although with an additional error.

The criteria for the stability of the periodic extension follow from the stability analysis of the associated DDE. For the nonlinear case the stability of the nonlinear DDE is equivalent to the stability of the linearized version of the original DDE. The basic result about the stability of a linear DDE is that any solution of a DDE is stable if and only if the largest real part of the roots of the characteristic equation is non-positive. The criteria for the roots of the characteristic equation to lie in the complex left half-plane have been studied extensively. A summary of the results can be found in .

## 5 Parameter Adaptation Dynamics

If the construction of periodic extensions is realized in the physiology of the brain then there must exist an algorithm that adjusts weights and/or delays to arrive at the needed values so that a population of neurons learns the appropriate weights after many repetitions of exposure to the external input. In this section we present examples of such learning algorithms. Of course we do not claim that either of the algorithm is actually implemented in the brain. Neither is learning a central point of this paper. Nevertheless it is instructive to consider some possibilities to estimate the difficulty of the problem.

We shall consider that the delays are known and fixed and only the weight matrices $`A_k`$ need to be determined. The delays can also be considered as variables instead or together with $`A_k`$ . Investigation of this possibility we leave for elsewhere.

A learning algorithm is not difficult to construct by introducing a Liapunov dynamics in the space of weights so that at the fixed point of the evolution we obtain the needed periodic extension. The weight dynamics can be consequently constructed by defining an error functional that measures the distance from the desired solution to the initial condition and using it as a Liapunov function. Consider, for example, ”slow” learning done during the multiple presentations - ”epochs” - of the same initial condition. The $`\mathrm{"}m\mathrm{"}`$th update of $`A_k`$ would then be done once per ”epoch”: one update for each time interval $`[m\tau _L,\left(m+1\right)\tau _L].`$ As an example we can take the $`\mathrm{"}m\mathrm{"}`$th epoch error functional as

$$V_m\left[A\right]=_{m\tau _L}^{\left(m+1\right)\tau _L}𝑑t\left|x\left(t\right)\varphi _0\left(t\right)\right|^2,$$
(41)
where $`x\left(t\right)`$ is the solution of the DDE with a given initial condition $`\varphi _0\left(t\right).`$ ”Fast” adaptation within the single epoch can also be defined in a similar way, although from the point of view of biology it might be more suitable for adaptation of delays. The error functional $`\left(\text{41}\right)`$ is phenomenologically appealing since both $`x\left(t\right)`$ and $`\varphi \left(t\right)`$ can be physically ”measured” by the encoding population of neurons. The discrete evolution in the space of weights with fixed delays can then be written as

$$\left(A_k\right)_{m+1}=\left(A_k\right)_{m+1}\epsilon \frac{}{A_k}V_m\left[A\right],$$
(42)
where $`\epsilon `$ is a small parameter. If the fundamental solution can be computed then $`x\left(t\right)`$ in $`\left(\text{41}\right)`$ is given by

$$x\left(t\right)=X\left(t\tau _L\right)\varphi \left(\tau _L\right)+\underset{k=1}{\overset{L}{}}_{\tau _L\tau _k}^{\tau _L}𝑑\tau X\left(t\tau \tau _k\right)A_k\varphi \left(\tau \right),$$

where the dependence of the fundamental solution $`X(t)`$ on the weights is given by (7,6). Computation of the appropriate derivatives of $`V[A]`$ involves computing the $`A`$derivatives of the fundamental solution. These can be written out as

$$\frac{X_{pq}}{\left(A_k\right)_{ij}}=\frac{1}{2\pi i}_{\left(c\right)}𝑑\lambda \mathrm{exp}\lambda \left(t\tau _k\right)\left(H^2\left(\lambda \right)\right)_{pi}\delta _{qj},$$
(43)
which enables us to compute $`\frac{}{A_k}V_m\left[A\right]`$.

A more direct approach to learning the weights is to use the characteristic equation itself. We know from the discussion in the preceding sections that matrices $`A_k`$ must be chosen in such a way that the characteristic equation has purely imaginary roots located at $`\lambda _n=\frac{2\pi i}{\tau _L}n,`$ $`n=Q,\mathrm{},Q.`$ Therefore an alternative Liapunov function for the weight dynamics can be chosen as

$`V_H\left[A\right]`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{n=Q}{\overset{Q}{}}}\left|detH\left(\lambda _n\right)\right|^2`$ (44)
$`H\left(\lambda \right)`$ $`=\lambda I{\displaystyle \underset{k=0}{\overset{L}{}}}A_k\mathrm{exp}\left(\lambda \tau _k\right)`$
The weight dynamics induced by this Liapunov function is given by

$$\dot{A}_k=\epsilon \mathrm{Re}\left(\underset{n=Q}{\overset{Q}{}}\mathrm{exp}\left(\lambda _n\tau _k\right)det\overline{H}\left(\lambda _n\right)H^T\left(\lambda _n\right)\right)$$
(45)
where bar denotes complex conjugation. This system of equations conceptually simpler than the one obtained from $`\left(\text{41}\right),`$ since it does not involve integration over time.

Yet another approach is to use the explicit form of the solutions for $`A_k`$ and define the norm in the space of matrices by

$`V_Y\left[A\right]`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{k=0}{\overset{L}{}}}tr\left(A_kY_k\right)\left(A_kY_k\right)^T`$ (46)
$`Y_k`$ $`={\displaystyle \underset{k=0}{\overset{2Q}{}}}T_{kn}^1\mathrm{\Gamma }_n`$
where $`T,\mathrm{\Gamma }`$ are defined by $`(\text{28}\text{29},\text{32},\text{33}).`$ Differentiation with respect to unknown $`A_k`$ results in

$$\dot{A}_k=\epsilon \left(\left(A_kY_k\right)\right)$$

The three algorithms we presented all apply only to the linear case. For applications in biology the DDE are typically nonlinear. In this case one can use a stochastic annealing approach to the error functional $`V_m\left[A\right]`$ defined by $`\left(\text{41}\right).`$ Another alternative is to restrict learning to the linear regime of the nonlinear DDE. This is possible because the typical nonlinearities in memory models, e.g. the Hopfield model, do have a region of definition of the nonlinearity where $`\sigma \left(z\right)z.`$

## 6 Summary

In this paper we presented a method for approximate storage of the trajectories of nonlinear dynamical systems as periodic or almost periodic solutions of the linear and nonlinear delay-differential equations.

Although the main motivation for this work was to develop a model for time-variable pattern storage by the brain, the results might have wider applicability. Indeed, they indicate that linear delay differential equations are, in a sense, universal models for nonlinear dynamical systems. Choosing appropriately large $`N`$ one in principle can store arbitrary (but finite) number of periodic orbits of a nonlinear dynamical system.

In addition to providing a method for storage of time-dependent patterns our model has other features that are attractive for a biological interpretation. For example, the method requires global synchronization of populations of neurons, since the initial conditions have to be represented via Fourier series with the same fundamental frequency $`\frac{2\pi }{\tau _L}`$. The maximum delay $`\tau _L`$ can be thought of the attention span of the population: periodic trajectories with period larger then $`\tau _L`$ cannot be stored. The model does not have unknown parameters. Once the delays are known the weights can be determined uniquely for a given trajectory of the environment.

There are a number of open questions about our method. Although the situation with the linear DDE seems to be clear, a more detailed investigation of the effect of the presence of nonlinearity should be carried out. One beneficial effect of nonlinearity could be multistability for a chosen set of the coefficients $`A_k,\tau _k`$: the existence of different, initial condition dependent periodic solutions similar to the existence of multiple fixed points in the Hopfield model. This would make multiple trajectory storage more efficient: a linear DDE can store additional orbits only by increasing its dimension linearly with the number of the orbits.

We have shown that nonlinearity brings along one interesting feature: multistability. More than one periodic solution can be stored for a given set of matrices $`A_k.`$ When nonlinearity is polynomial, one can estimate that maximum possible number of memories grows exponentially in the number of Fourier coefficients. Additional indication that the multistability should exist comes from considering a continuum limit of $`\left(\text{35}\right),`$ a limit that is reasonable to take when there are large number of delays that depend smoothly on their index. This is the situation when neurons in the brain have large dendritic trees. With obvious definitions the continuum limit can be formally written as

$$\dot{\varphi }\left(t\right)=\lambda _0^1𝑑sA\left(s\right)\sigma \left(\varphi \left(t\tau \left(s\right)\right)\right)$$
(47)
for $`\lambda `$ a parameter and for some monotonic delay function $`\tau \left(s\right).`$ Here we normalized the matrix-valued function $`A\left(s\right)`$ so that $`tr\left(A\left(s\right)^TA\left(s\right)\right)=1.`$ With the use of the Fourier transform this equation can be related to the Volterra and Hammerstein classes of nonlinear integral equations, which exhibit bifurcations in the parameter $`\lambda `$. The theory of nonlinear integral equations has been exhaustively studied and a number of criteria are available for determining the bifurcation points . These depend mainly on the growth properties of the nonlinearity $`\sigma \left(z\right)`$. A more detailed analysis of $`\left(\text{47}\right)`$ shall be considered elsewhere.

###### Acknowledgement 1

I would like to thank Holger Kantz and members of the Nonlinear Time Series Analysis Group at MPI-PKS for useful discussions.


format: text

text:
# Dynamically tuning away the cosmological constant in effective scalar tensor theories

## 1. Introduction:

It is now widely accepted that a theory of gravitation described by the Einstein - Hilbert action:

$$S=d^4x\sqrt{g}[\frac{1}{8\pi G}R+\mathrm{\Lambda }]$$
(1)
fares wonderfully well to classical precision tests. The success is impressive enough to rule out most alternative theories at the classical level . However, the theory described by eqn(1) has severe theoretical inconsistencies. The appearance of a dimensional gravitational coupling constant $`G`$, and the smallness of the cosmological constant $`\mathrm{\Lambda }`$, are two pathologies that impede any attempt to treat this action as a fundamental quantum theory. The dimensional coupling is primarily responsible for the non - renormalizability of the theory  while the smallness of $`\mathrm{\Lambda }`$ can not be justified in any reasonable theoretical model  Besides these theoretical problems, a cosmological model based on the above action has to contend with dynamical inconsistencies: reasonable initial conditions can not be dynamically realized on account of the horizon and fine - tuning (flatness) problems. A resolution of these dynamical inconsistencies using inflation, without addressing the cosmological constant problem is, at best, incomplete. Inflation is founded upon the possibility that the early universe witnessed several phase transitions - each producing a large change in the effective value of the cosmological constant. Be it so, the last transition is required to exit to a state of a tiny cosmological constant. All versions of inflation require fine tuning of model parameters to get inflation going and / or to get inflation exit gracefully to the standard early hot universe . The least troublesome version of inflation, the chaotic model, assigns to $`\mathrm{\Lambda }`$ a vanishing value by hand. Anthropic arguments are invoked to justify a small value for $`\mathrm{\Lambda }`$  and have invited expressions of despondency and a need for a better ansatz. We share Witten’s  reservations on the use of anthropic arguments in general.

Over the last two decades there is an emerging consensus that the Einstein - Hilbert action ought to be regarded, at most, as a correct low energy limit of a decent quantum theory. Superstring theory has emerged as a promising candidate in this respect. It is hoped that an effective gravity theory would emerge from the dynamics of an appropriate string theory. A lack of understanding of the vanishing (or the tininess) of the cosmological constant is the key obstacle to making realistic string-inspired particle physics models. A convincing mechanism of supersymmetry breaking has not emerged as yet. All known approaches generate a large cosmological constant and must therefore be regarded as incorrect .

On the other hand, general considerations regarding the structure of string theory suggest that general relativity may acquire modifications, even at energies considerably lower than the Planck scale. Dilaton gravity theories in ten dimensions have recently played a prominent role in a search for a consistent effective low energy limit emerging from string inspired quantum theory of gravity . This has resulted in a resurrection of interest in scalar \- tensor and Brans-Dicke models. Such models have a characteristic coupling of the ricci scalar curvature with a function of a scalar field - defining a non - minimal coupling (NMC). A search for a suitable compactification scheme that would yield the Einstein - Hilbert action in four dimensions is on within the framework of such models. To achieve this, the foremost requirement is to have a provision for an appropriate fixing of NMC to give a universal gravitational constant.

The NMC could, for example, be determined by the minimum of an effective potential of the scalar field. Such a potential may arise as a result of the compactification scheme itself or as a result of quantum corrections . Determining the NMC by the minimum of an effective potential leads to a realization of a theory indistinguishable from the Einstein - Hilbert theory at low energies . However, such a scheme falls short of addressing itself to the cosmological constant problem as again, there is nothing to ensure the tininess of the minimum value of the effective potential.

There have been attempts to dynamically tune the effective cosmological constant to a vanishingly small value in classes of NMC theories . All such attempts are frustrated by a “no-go theorem” . In these schemes, the scalar field develops an instability, the non-minimal coupling diverges and the scalar field stress tensor evolves to neutralize the effective value of the cosmological constant to zero. The trouble with such a scheme is that with the diverging NMC, the effective gravitational constant itself vanishes. This no-go situation has not been successfully evaded. One has not found a natural way to stabilize the NMC to a bounded value that could be identified with $`(8\pi G)^1`$ and, at the same time, tune the cosmological constant to a small value.

In a sweeping review, Weinberg  has made out a case for the use of a judicious anthropic argument as an explanation for the tininess of the cosmological constant. One considers formalisms that allow an effective cosmological constant to appear dynamically in the action and appeals to some unique feature that would pick out a small value for this constant in the functional integral. The most explored principal is the imposition of a conditional probability on the functional integral that would allow for a large enough universe to exist and also allow for gravitational clustering at a sufficiently small redshift in the universe. This imposes severe constraints on the cosmological constant. Indeed the most recent attempts  to justify a small $`\mathrm{\Lambda }`$ uses an anthropic argument to have a vanishing minimum to an effective “tracker” potential for a scalar field. The peculiar profile of such potentials, and a slow asymptotic rolling of the scalar field near the minimum, can in principle allow for an effective $`\mathrm{\Lambda }`$ dominated cosmology at the current epoch.

A feature that would pick out a small cosmological constant using an objective physical principle rather than having to take a recourse to an anthropic principle, is worth exploring. In this article we describe solutions in a class of non-trivial NMC theory having a minimum energy at a vanishing cosmological constant. The solutions arose in our search for an ansatz that would allow for a NMC with a non-trivial spatial dependence. As a result of a scalar field coupling to a fermion field, the expectation value of the scalar field can acquire a non-trivial value inside compact domains. But for the non-minimal coupling (of the scalar field) to the scalar curvature, these configurations are similar in character to the non-topological solutions suggested by Lee and Wick . With NMC, one gets a peculiar spatial variation of effective gravity. We shall require $`G`$ to be held to a non-trivial value inside compact domains. This is sufficient to have long distance gravitational effects indistinguishable from those expected from general theory of relativity (GTR). We point out two ways this could be possible: (a) the domains could be large enough to accommodate structures as large as a typical galaxy or (b) the compact domains could be of microscopic (sub-nuclear) size. This is described in the next section where we also describe some essential aspects of a generalized effective scalar tensor theory. The action is chosen to include scalar functions multiplying the ricci scalar as well as a fermionic part. Stability of large fermion number non-trivial solutions requires their interior to be near the zero of the effective potential.

In the last section we discuss features of a toy cosmological model that such a theory can support. The struture has definitive, falsifiable predictions. In standard cosmology, one has pinned down hopes on special initial conditions, motivated by anthropic considerations to account for the observed cosmological constant. Such tautologies, when applied to the universe as a whole, have little predictive power. In the structure described in this article, any fine tuning of coupling parameters of the fermion lagrangian that may be required to have stable gravity domains of required size, has testable consequences.

## 2. A generalized scalar tensor theory:

### 2(a). The action and field equations:

Consider a classical theory of a scalar field $`\varphi `$ coupled to the scalar curvature $`R`$, through an arbitrary function $`U(\varphi )`$, in the action:

$$S=d^4x\sqrt{g}[U(\varphi )R+\frac{1}{2}^\mu \varphi _\mu \varphi V(\varphi )+L_m]$$
(2)
We shall treat $`S`$ as a classical, effective, phenomenological action that could be good enough to do tree level (“on-shell”) physics. If need arises, we shall consider non polynomial forms for the effective potential and the NMC: $`V(\varphi )`$ and $`U(\varphi )`$. Such profiles have been used as “tracker-field” potentials  and can be justified in particle physics models with dynamic symmetry breaking or non-perturbative effects. In particular they can arise from non-perturbative effects that lift flat directions in supersymmetric gauge theories as well as in modulii fields in string theory  (see also ). Strictly speaking, we would consider it premature to justify the forms for these functions at this formative stage on the basis of fundamental physics. For our purpose, we simply regard eqn(2) as the action of a parametrized Brans-Dicke fluid. Any arbitrary (cosmological) constant in the theory is included in $`V(\varphi )`$. $`L_mL_w+L_\psi `$ is the contribution from the rest of the matter fields. It includes a contribution from a Dirac fermion field:

$$L_\psi \overline{U}(\varphi )^1[\frac{1}{2}[\overline{\psi }\stackrel{}{D}_\mu \gamma ^\mu \psi \overline{\psi }\gamma ^\mu \stackrel{}{D}_\mu \psi ]m(\varphi )\overline{\psi }\psi ]$$
(3)
Here $`D_\mu `$ is the spin covariant derivative . This Fermion action has an overall scaling function $`\overline{U}(\varphi )^1`$. Such an arbitrary scaling can be absorbed in a redefinition of the Dirac field: $`\psi \psi ^{}\psi /\overline{U}^{\frac{1}{2}}`$ provided $`\overline{U}`$ were bounded. Indeed, by such a scaling (for a $`\mathrm{𝑏𝑜𝑢𝑛𝑑𝑒𝑑}\overline{U}`$), the $`\varphi `$ dependence in $`L_\psi `$ can be constrained to the mass munction $`m(\varphi )`$. The anti-symmetrized derivative appearing in the fermion lagrangian ensures the cancellation of the derivative of a real, non-vanishing, bounded scaling function. 111This result is not specific to the chosen form of the Dirac Lagrangian. Covariant spinor fields of arbitrary weights can be defined as geometrical objects over (direct products of) two dimensional complex spaces. An arbitrary, non vanishing, real scaling of such a spinor can be absorbed into a real part of the trace of a spinor connection and does not contribute to the covariant derivative of the spinor and therefore has no interaction with other fields . If $`\overline{U}(\varphi )^1`$ vanishes at any point, the Dirac action identically vanishes. There is no dynamics for the Dirac field at such a point. The invariance of the Dirac Lagrangian under arbitrary phase change: $`\delta \psi =iϵ\psi `$ implies the vanishing divergence of the current:

$$J^\mu =\overline{U}(\varphi )^1\overline{\psi }\gamma ^\mu \psi $$
$`(4)`$
and the conservation of the charge:

$$Q_\mathrm{\Sigma }\overline{U}(\varphi )^1\overline{\psi }\gamma ^o\psi $$
$`(5)`$
The Dirac particles would be confined to regions where $`\overline{U}`$ is bounded i.e. $`\overline{U}(\varphi )^1`$ is non - vanishing. The only effect of a scaling function $`\overline{U}(\varphi )^1`$ that vanishes outside a compact domain is to confine the fermions within the domain. Inside the domain the scaling can be identically absorbed in a rescaling of the fermion field. In what follows we shall essentially do this. This leads to conditions similar to those that arise in non-topological soliton solutions in the Lee-Wick model . One could just as well have followed the Lee-Wick model in which confinement occurs for fermions that are not on-shell outside a domain on account of a large coupling with the Higgs field.

Our foremost task would be to demonstrate the vanishing of the covariant divergence of the stress tensor of the rest of the matter fields described by $`L_w`$ in eqn(2). Requiring the action to be stationary under variations of the metric tensor and the fields $`\varphi \&\psi `$ gives the equations of motion:

$$U(\varphi )[R^{\mu \nu }\frac{1}{2}g^{\mu \nu }R]=\frac{1}{2}[T_w^{\mu \nu }+T_\varphi ^{\mu \nu }+\mathrm{\Theta }^{\mu \nu }+2U(\varphi )^{;\mu ;\nu }2g^{\mu \nu }U(\varphi )]_{;\lambda }^{;\lambda }]$$
$`(6)`$
$$g^{\mu \nu }\varphi _{;\mu ;\nu }+\frac{V}{\varphi }R\frac{U}{\varphi }+\frac{m}{\varphi }\overline{\psi }^{}\psi ^{}=0$$
$`(7)`$
$$\gamma ^\mu D_\mu \psi ^{}+m(\varphi )\psi ^{}=0$$
$`(8)`$
$$D_\mu \overline{\psi }^{}\gamma ^\mu m(\varphi )\overline{\psi }^{}=0$$
$`(9)`$
Here $`T_w^{\mu \nu }`$, $`\mathrm{\Theta }_\nu ^\mu `$ are the energy momentum tensors constructed from $`L_w`$ and $`L_\psi `$ respectively, and

$$T_\varphi ^{\mu \nu }^\mu \varphi ^\nu \varphi g^{\mu \nu }[\frac{1}{2}^\lambda \varphi _\lambda \varphi V(\varphi )]$$
$`(10)`$
$`L_w`$ is taken to be independent of $`\varphi `$ and $`\psi `$. The two Dirac equations (8) and (9) ensure $`L_\psi `$ to be null. The scalar field equations are therefore independent of $`\overline{U}(\varphi )`$. The Fermion stress tensor is simply:

$$\mathrm{\Theta }_\nu ^\mu \frac{1}{2}[\overline{\psi }^{}\stackrel{}{D}_\nu \gamma ^\mu \psi ^{}\overline{\psi }^{}\gamma ^\mu \stackrel{}{D}_\nu \psi ^{}]$$
$`(11)`$
The covariant divergence of (6) is easily seen to reduce to:

$$\mathrm{\Theta }_{\nu ;\mu }^\mu =\frac{m}{\varphi }_\nu \varphi \overline{\psi }^{}\psi ^{}$$
$`(12)`$
Thus there is a violation of equivalence principle as far as the Dirac field is concerned. However, in a region where the scalar field gradient vanishes or for a $`\varphi `$ independent fermion mass, the covariant divergence of the fermion field stress tensor vanishes. To see how the equivalence principle strictly holds for the rest of the matter fields, i.e.: $`T_{w;\nu }^{\mu \nu }=0`$, consider the covariant divergence of eqn(6). From the contracted Bianchi identity satisfied by the Einstein tensor, we obtain

$$U(\varphi )_{,\nu }[R^{\mu \nu }\frac{1}{2}g^{\mu \nu }R]=\frac{1}{2}[T_{w;\nu }^{\mu \nu }+t_{;\nu }^{\mu \nu }+\mathrm{\Theta }_{;\nu }^{\mu \nu }]$$
$`(13)`$
with

$$t^{\mu \nu }T_\varphi ^{\mu \nu }+2U(\varphi )^{;\mu ;\nu }2g^{\mu \nu }U(\varphi )_{;\lambda }^{;\lambda }$$
$`(14)`$
Using the identity:

$$U(\varphi )^{;\rho }R_{\rho \alpha }=U(\varphi )_{;\lambda ;\alpha }^{;\lambda }U(\varphi )_{;\alpha ;\lambda }^{;\lambda }$$

and eqn(12), eqn(13) reduces to

$$\frac{1}{2}U(\varphi )^{,\mu }R=\frac{1}{2}[T_{w;\nu }^{\mu \nu }+T_{\varphi ;\nu }^{\mu \nu }+^\mu m\overline{\psi }^{}\psi ^{}]$$
$`(15)`$
Finally, using the equation of motion for the scalar field (7), all the $`\varphi `$ dependent terms cancel the left hand side - giving the vanishing of the covariant divergence of the (w-) matter stress energy tensor. This in turn ensures that equivalence principle holds \- in turn assuring the geodesic law.

### 2(b). The conserved energy:

Next, we write down the expression for a conserved pseudo energy momentum tensor . Defining

$$\mathrm{A}\sqrt{g}g^{\sigma \rho }[\mathrm{\Gamma }_{\sigma \rho }^\alpha \mathrm{\Gamma }_{\alpha \beta }^\beta \mathrm{\Gamma }_{\beta \rho }^\alpha \mathrm{\Gamma }_{\alpha \sigma }^\beta ]$$
$`(16)`$
$$\mathrm{B}[U\mathrm{A}\sqrt{g}g^{\sigma \rho }\mathrm{\Gamma }_{\sigma \alpha }^\alpha U_{,\rho }+\sqrt{g}g^{\sigma \rho }\mathrm{\Gamma }_{\sigma \rho }^\alpha U_{,\alpha }]$$
$`(17)`$
and $`\widehat{\mathrm{B}}\mathrm{B}+\sqrt{g}L_{\varphi +\psi }`$, the conserved pseudo energy momentum vector defined over a hypersurface $`\mathrm{\Sigma }`$ is:

$$P_\mu _\mathrm{\Sigma }𝑑\mathrm{\Sigma }[\sqrt{g}T_{w\mu }^o\widehat{\mathrm{B}}\delta _\mu ^o\frac{\widehat{\mathrm{B}}}{g_{,o}^{\tau \beta }}g_{,\mu }^{\tau \beta }\frac{\widehat{\mathrm{B}}}{\varphi _{,o}}\varphi _{,\mu }\overline{\psi }_{,\mu }\frac{\widehat{\mathrm{B}}}{\overline{\psi }_{,o}}\frac{\widehat{\mathrm{B}}}{\psi _{,o}}\psi _{,\mu }]$$
$`(18)`$
This is also expressible as a surface integral over a 2 dimensional boundary of $`\mathrm{\Sigma }`$:

$$P_\mu =(\frac{\widehat{\mathrm{B}}}{g_{,j}^{\mu \nu }}g^{0\nu })𝑑\mathrm{\Sigma }_j$$
$`(19)`$
Thus in a generalised Brans-Dicke theory, the generalised energy momentum is determined by the metric tensor and its derivatives on a 2-dimensional boundary.

### 2(c). The dynamic equations:

We consider classical solutions to eqns.(6-9) for a fixed number of fermions. As demonstrated in  and  the thermodynamic potentials for the fermion field are assumed to be adequately described in terms of the chemical potential and the temperature. Denoting the fermion density by $`S_f`$, eqn(7) reads:

$$\mathrm{}\varphi +\frac{V}{\varphi }R\frac{U}{\varphi }+\frac{m}{\varphi }S_f=0$$
$`(20)`$
The trace of eqn(6) gives:

$$U(\varphi )R=(3U^{\prime \prime }(\varphi )+\frac{1}{2})\varphi ^{,\alpha }\varphi _{,\alpha }+2V(\varphi )+\frac{1}{2}mS_f3U^{}(\varphi )\mathrm{}\varphi $$
$`(21)`$
Substituting in eqn(20) gives:

$$\mathrm{}\varphi +\frac{U^{}(\varphi )(3U^{\prime \prime }(\varphi )+\frac{1}{2})}{1+3\frac{U^2}{U}}\varphi ^{,\alpha }\varphi _{,\alpha }+\frac{V^{}+m^{}S\frac{U^{}}{U}(2V+\frac{1}{2}mS)}{1+3\frac{U^2}{U}}$$
$`(22)`$
### 2(d). An example of tuning of $`\mathrm{\Lambda }`$:

Consider a familiar NMC theory with $`U(\varphi )=\beta /8\pi \xi \varphi ^2/2`$ and $`V(\varphi )\mathrm{\Lambda }_0/8\pi `$ a constant. We define natural units in which $`\beta `$ takes on a unit value. The field equations simply read:

$$\mathrm{}\varphi +\xi R\varphi =0$$
$`(23)`$
and

$$G_{\mu \nu }+\mathrm{\Lambda }_Og_{\mu \nu }=8\pi T_{\mu \nu }$$
$`(24)`$
with

$$T_{\mu \nu }=\varphi _{,\mu }\varphi _{,\nu }1/2g_{\mu \nu }\varphi _{,\rho }\varphi ^{,\rho }\xi \varphi ^2(R_{\mu \nu }1/2g_{\mu \nu }R)+\xi g_{\mu \nu }\mathrm{}\varphi ^2\xi \varphi _{}^{2}{}_{;\mu ;\nu }{}^{}$$
$`(25)`$
In terms of $`\mathrm{\Phi }\varphi ^2`$, the scalar field eqn(23) reads

$$\mathrm{}\mathrm{\Phi }=2\xi R\mathrm{\Phi }+\frac{1}{2\mathrm{\Phi }}_\mu \mathrm{\Phi }^\mu \mathrm{\Phi }$$
$`(26)`$
while the trace of eq(24) is simply:

$$R+4\mathrm{\Lambda }_0=4\pi [\frac{\mathrm{\Phi }_{,\rho }\mathrm{\Phi }^{,\rho }}{2\mathrm{\Phi }}+2\mathrm{\Phi }R\xi +6\xi \mathrm{}\mathrm{\Phi }]$$
$`(27)`$
The wave equation for the scalar field thus reduces to:

$$\mathrm{}\mathrm{\Phi }=\frac{4\mathrm{\Lambda }_0R}{4\pi (16\xi )}$$
$`(28)`$
We shall follow the dynamics of $`\varphi `$ in a FRW metric:

$$ds^2=dt^2a^2(t)[\frac{dr^2}{1kr^2}+r^2(d\theta ^2+sin^2\theta d\phi ^2)]$$
$`(29)`$
For a homogeneous, time dependent scalar field, the above equations reduce to:

$$\frac{3}{8\pi }[(\frac{\dot{a}}{a})^2+\frac{k}{a^2}]\frac{\mathrm{\Lambda }_0}{8\pi }=\frac{\dot{\mathrm{\Phi }}^2}{8\mathrm{\Phi }}+3\xi [(\frac{\dot{a}}{a})^2+\frac{k}{a^2}]\mathrm{\Phi }+3\xi (\frac{\dot{a}}{a})\dot{\mathrm{\Phi }}$$
$`(30)`$
and

$$\ddot{\mathrm{\Phi }}+3(\dot{a}/a)\dot{\mathrm{\Phi }}=\frac{4\mathrm{\Lambda }_0R}{4\pi (16\xi )}$$
$`(31)`$
We look for solutions with the scale factor having a power law evolution: $`a(t)=ht^\alpha `$, with $`h,\alpha `$ arbitrary constants. The equation for $`\mathrm{\Phi }`$ integrates exactly for all $`\alpha 1/3`$ in terms of integration constants $`K_o,C`$:

$$\mathrm{\Phi }=At^2/2+K_0t^{13\alpha }+C$$
$`(32)`$
where

$$A=\frac{\mathrm{\Lambda }_0}{\pi (16\xi )(1+3\alpha )}$$
$`(33)`$
(For $`\alpha =1/3`$ one has $`K_0ln(t)`$ instead of the second term. From what just follows, this is not feasible as a late time solution). On substituting in eqn(20) it is found that it is possible to have $`k0`$ solutions for large times only if $`\alpha 1`$. For $`\alpha >1`$, the solutions are independent of the curvature constant (i.e. hold for all $`k`$). The scale $`\alpha `$ directly relates to the parameter $`\xi `$ as: $`\xi =(4\alpha 2)^1`$. Similarly, for $`\alpha =1`$, $`\xi =(2+2k/h)^1`$. For all these cases, scalar condensates develop to have their energy cancell the effect of $`\mathrm{\Lambda }_0`$:

$$T_{\mu \nu }\frac{1}{8\pi }\mathrm{\Lambda }_0g_{\mu \nu }+O(t^2)$$
$`(34)`$
The cancellation occurs over time scales large as compared to the time scales determined by dimensional parameters in the action. If one merely uses typical particle physics time scales, this implies that the stress energy of the scalar field would compensate the cosmological constant for the entire history of the universe except the very early history. It is the coupling of the scalar field to the spacetime curvature that is primarily responsible for the above instability of the scalar field which in turn leads to an efficient damping of the cosmological constant. Unfortunately, the effective value of the gravitational constant goes as $`G_{eff}=(\beta +8\pi |\xi |\varphi ^2)^1`$. This becomes unacceptably small at large times. This is an example of a “no-go” situation . This behaviour was noted in  for a flat $`k=0`$ model. We conclude that the presence of spatial curvature does not alter the results. It is also straightforward to demonstrate that the results are independent of the presence of the Planck sacle determined by $`\beta ^18\pi G_0M_{plank}^2`$. Indeed in the absence of $`\beta `$ in the NMC function $`U(\varphi )`$, the FRW scale evolves as above: $`a(t)t^\alpha `$ with $`\alpha `$ constrained as before in the presence of $`\beta `$, the stress tensor rapidly cancells the cosmological constant while the effective gravitational constant again approaches zero over time scales large in comparison to the particle physics time scales.

The above behaviour of a power law behaviour of the FRW scale factor and the cancellation of the vacuum energy has been generalized to scalar tensor theories within a large class of scalar field effective potentials . There have also been attempts to use the above results to constrain dimensional parameters in the theory (see eg ). One equates the observed gravitational constant to that calcuted above using conservative age estimates of the universe. Early universe nucleosynthesis constraints and recent Viking radar experiments severly constrain the NMC in a scalar tensor theory. The particle physics parameters that follow turn out to be rather unnatural. Instead we propose and explore the possibility that the above analysis successfully accounts for a NMC diverging with time and that an effective gravitational constant over most of the universe is indded vanishing. However, this is not the constant which ought to be equated with the Newtonian gravitational constant. Suspecting the problem to lie in the assumption of homogeneity, it has been suggested  to look for a stabilization of $`G_{eff}`$ in an inhomogeneous model.

### 2(e). A strategy for effective gravity:

Taking a cue from the above, we consider a theory described by eqns(2 & 3). We shall be interested in dynamic NMC’s that diverge, as described above, outside bounded domains. It is convenient to recast the equations in terms of the effective planck length parameter which goes to zero as the NMC diverges. One such reparametrization of the scalar field is $`\varphi M^2/\varphi `$. The NMC becomes $`U(\varphi )\beta \xi M^4/\varphi ^2`$ . With this reparametrization, the results of the last section translate to the vanishing of the effective gravitational and cosmological constants together with the vanishing of $`\varphi `$ over time scales large in comparison to any particle physics time scales in the theory. In other words we consider an equivalent problem that has the NMC divergent at any point that we take as $`\varphi =0`$ without any loss of generality. Thus $`U(\varphi =0)\mathrm{}`$. The model has a flat, Minkowski spacetime as a fixed point global solution with $`\varphi =0`$ for an arbitrary equation of state of matter. However, there are non-trivial solutions to the model as well . We are particularly interested in solutions that have $`\varphi 0`$ outside compact domains but not in the interior. This would mean that the gravitational constant would vanish in the exterior. However, this does not imply doing away with gravitation between such domains. To see this, consider the Einstein tensor that follows from equations (2) and (6):

$$G_{\mu \nu }=[U(\varphi )]^1\times sourceterms$$
$`(35)`$
From the expansion of the metric about a flat space metric: $`g_{\mu \nu }=\eta _{\mu \nu }+h_{\mu \nu }`$, in appropriate (harmonic) coordinates, one gets:

$$\mathrm{}h_{\mu \nu }=[U(\varphi )]^1\times sourceterms$$
$`(36)`$
The retarded solution reads:

$$h_{\mu \nu }(x,t)=\frac{d^3x^{}}{|xx^{}|}[U(\varphi )]^1\times sourceterms(x^{},t|xx^{}|)$$
$`(37)`$
Thus as long as the effective gravitational constant at the retarded source point $`(x^{},t|xx^{}|)`$ is bounded, the vanishing of the effective gravitational constant at the field point is of no consequence for the above solution. All one has to discover is the means of having the NMC locked to a universal value in non-trivial domains. The mechanism we invoke is to have the non-minimally coupled scalar field locked at the minimum of its effective potential .

### 2(f). Large Non trivial solutions:

We recall studies on non-trivial configuration that have $`\varphi `$ vanishing outside, and $`N_f>>0`$ $`\psi `$ fermions trapped inside, compact domains in the Lee-Wick model . With the fermion number conserving, the minimum energy configuration of such solutions are non-topological soliton solutions \[NTS’s\] in the theory. We recite the essential properties of such an NTS. Lee and Wick considered a theory with a Lagrangian:

$$L=\frac{1}{2}_\mu \varphi ^\mu \varphi V(\varphi )+\overline{\psi }\gamma ^\mu _\mu \psi m\overline{\psi }[1\frac{\varphi }{\varphi _0}]\psi $$
$`(38)`$
where

$$V(\varphi )=\frac{1}{2}m_\varphi ^2\varphi ^2[1\frac{\varphi }{\varphi _0}]^2$$
$`(39)`$
With Fermion number conserved in this theory, solutions representing a degenerate distribution of fermions with a total fermion number $`N_f`$, at $`\varphi =\varphi _0`$, would remain confined there if the total energy of a fermion is less than its on-shell energy at $`\varphi =0`$. The NTS is a ball of radius $`r_b`$ inside which the scalar field is locked at $`\varphi \varphi _0`$. Across the surface of the ball of thickness $`m_\varphi ^1`$, the scalar field transits to $`\varphi =0`$. As a matter of fact, the spherical shape of the distribution follows from degeneracy of the fermion ensemble. For a degenerate fermion distribution, the fermion density is determined in terms of the chemical potential that occurs as a lagrange multiplier for the conserved fermi number of the ensemble. This determines the volume for a fixed conserved number of degenerate fermions. A positive definite surface energy would break the degeneracy of volume preserving (i.e. fermion number conserving) deformations of the distribution \- choosing a configuration of minimum surface area. This gives the minimum energy configuration to be a sphere. In terms of the chemical potential $`\mu `$, the fermion density goes as $`\rho _f\mu ^4`$. The fermion number density goes as $`n_f\mu ^3`$. With the total energy and fermion number: $`E_f=\rho _fv\mu ^4r_b^3`$, $`N_f=n_fv\mu ^3r_b^3`$, the chemical potential eliminates to give $`E_fN_f^{4/3}/r_b`$. For a large enough $`r_b`$, the surface tension is given by the expression: $`sm_\varphi \varphi _0^2/6`$, the total surface energy being $`E_s=4\pi sr_b^2`$. For a given $`N_f`$, one may vary the radius to seek the minimum of the volume and surface energies to give $`E_f=2E_s`$. The mass of the soliton is expressed as $`M=E_f+E_s=3E_s=3E_f/2=12\pi sr_b^2`$. This allows one to express:

$$N_fs^{3/4}r_b^{9/4};Ms^{1/3}N_f^{8/9};r_bs^{1/3}N_f^{4/9}$$
$`(40)`$
Thus with a given surface tension, one could get arbitrarily massive NTS’s by increasing the fermion number. However, the above flat spacetime analysis can not be extended past the Schwarzschild limit. In the Lee-Wick model, the Schwarzschild limit can be determined as the criteria for stability against gravitational collapse. As a matter of fact, gravitational instability occurs even before this limit is reached . We shall use the Schwarzschild limit as a benchmark to keep sufficiently clear-off and as a convenience that ensures that the flat spacetime analysis holds. This limit is simply obtained by equating the radius of the NTS to the Schwarzschild radius for a critical mass ball $`M_c`$: equating $`r_b=2GM_c`$, and with the expression for the mass of a ball, one gets $`M_c(48\pi G^2s)^1`$. For $`s(30Gev)^3`$, $`M_c`$ is roughly $`10^{15}M_{}`$ with the critical radius some $`10^2`$ Pc.

One can have NTS’s as large as some 10’s of kilo parsecs. This would require a class of low mass fermions that could keep the NTS from collapsing. In a hot big bang cosmology, any species of particles that decoupled from equilibrium when it was relativistic, ought to have the same relic density as say the relic photons or neutrinos: some $`200`$ per cc. To get $`10`$ to $`100`$ KPc NTS’s, it would suffice to have $`s(Mev)^3`$. This gives $`N_f10^{72}`$ to $`10^{75}`$ with the mass of the ball some 12 orders of magnitude smaller than the the Schwarzschild (critical) mass for $`s(Mev)^3`$. We conclude that for such a surface tension and for NTS’s even as large as hundreds of kilo parsecs, one can consistently ignore the spacetime curvature in such solutions. In the following, we do this whenever convenience demands.

In the model described by eq(3), fermions $`\psi `$ are confined to $`\varphi 0`$ domains. It is easy to follow  to demonstrate the existence of non-trivial solutions. In the “weak (gravitational) field approximation”, that would justify retaining only a first order deviation from a flat metric, the metric can be expressed in terms of the spherical \[Schwarzschild\] coordinates:

$$ds^2=e^{2u}dt^2e^{2\overline{v}}dr^2r^2[d\theta ^2+sin^2\theta d\phi ^2]$$
$`(41)`$
We look for a static solution describing the scalar field trapped to a value $`\varphi =\varphi _{in}`$ in the interior of a sphere of radius $`r_b`$ and making a transition across a thin surface to $`\varphi =0`$ outside. The fermi gas trapped inside the soliton is described, in the Thomas Fermi approximation , by the following distribution in momentum space: $`n_k=\theta (kk_f)`$, $`k`$ being the momentum measured in an appropriate local frame that depends on $`r`$, and $`k_f`$ the fermi momentum. The fermion energy density is given by:

$$W=\frac{2}{8\pi ^3}d^3kn_kϵ_k$$
$`(42)`$
with $`ϵ_k=\sqrt{k^2+m(\varphi _{in})^2}`$. The fermion number density $`\nu _f`$ and the non-vanishing components of fermion stress energy tensor are:

$$\nu _f=\frac{2}{8\pi ^3}d^3kn_k$$
$`(43)`$
$$T_t^t=W$$

$$T_r^r=T_\theta ^\theta =T_\phi ^\phi =T_\rho ^\rho T=\frac{2}{8\pi ^3}d^3kn_k\frac{k^2}{3ϵ_k}$$
$`(44)`$
The trace of the stress tensor is just:

$$T_\mu ^\mu =W3T=m(\varphi _{in})S$$
$`(45)`$
with S the scalar density:

$$S=\frac{2}{8\pi ^3}d^3k\frac{n_k}{ϵ_k}m(\varphi _{in})$$
$`(46)`$
Defining $`8\pi G_{in}U(\varphi _{in})^1`$ as the effective interior “gravitational constant”, the metric field equation in the interior can be expressed in the above spherical coordinates as:

$$r^2G_t^t=e^{2\overline{v}}12e^{2\overline{v}}r\frac{d\overline{v}}{dr}=8\pi G_{in}r^2[W+V(in)]$$
$`(47)`$
$$r^2G_r^r=e^{2\overline{v}}1+2e^{2\overline{v}}r\frac{du}{dr}=8\pi G_{in}r^2[TV(in)]$$
$`(48)`$
$$r^2G_\theta ^\theta =e^{2\overline{v}}[r^2\frac{d^2u}{dr^2}+[1+r\frac{du}{dr}]r\frac{d}{dr}(u\overline{v})]=8\pi G_{in}r^2[TV(in)]$$
$`(49)`$
The scalar field satisfies:

$$\varphi _{;\mu }^{;\mu }+V^{}+m^{}(\varphi )SU^{}R=0$$
$`(50)`$
Taking the trace of the Einstein tensor $`G_\mu ^\mu `$, the Ricci scalar $`R`$ substitutes in eqn(24) to give the following radial equation for the scalar field in the weak field limit:

$$\frac{d^2\varphi }{dr^2}+\frac{2}{r}\frac{d\varphi }{dr}+F(U)[\frac{d\varphi }{dr}]^2=\frac{dW}{d\varphi }$$
$`(51)`$
where

$$\frac{dW}{d\varphi }[V^{}+m^{}(\varphi )S\frac{U^{}}{U}(\frac{mS}{2}+2V)]/[1+\frac{3U^2}{U}]$$
$`(52)`$
and

$$F(U)\frac{U^{}}{U}(\frac{1}{2}+3U^{\prime \prime })/[1+\frac{3U^2}{U}]$$
$`(53)`$
A sufficient condition for existence for a large NTS with $`\varphi _{in}`$ constant is the vanishing of $`W^{}`$ at $`\varphi _{in}`$. A convenient way to arrange this is to have $`V(\varphi )`$ to have a minima at $`\varphi _{in}`$, where $`U^{}/U`$ is small, and to have fermions massless at $`\varphi _{in}`$. The zero of $`W^{}`$ then coincides with that for $`V^{}`$. For a profile of $`V`$ described in Fig(1), and $`U(\varphi )=M_p^2+M^4/\varphi ^2`$ (with $`\varphi ^2>>M^4/M_p^2`$), the profile of $`W`$ closely follows that of $`V`$ except near $`\varphi 0`$. The divergence of the NMC at $`\varphi =0`$ ensures the vanishing of $`W^{}(\varphi =0)`$. The divergence of the NMC and in turn that of $`F(U)`$ also ensures that $`\varphi =0`$ is a solution with the Ricci scalar vanishing for arbitrary $`V(\varphi =0)`$. A typical non-trivial NTS solution has the scalar field locked to $`\varphi =0`$ outside a compact spherical domain of radius $`r_b`$ and transits to $`\varphi =\varphi _{in}`$ across the boundary:

$$1\frac{\varphi }{\varphi _{in}}exp[V^{\prime \prime }(\varphi _{in})(r_br)]$$

The transition zone is thus related only to $`V^{\prime \prime }(\varphi _{in})`$ for large $`r_b`$. The radius $`r_b`$ is determined by the energetics of the ball and is governed by the number of trapped fermions inside the ball. The pressure of the degenerate fermions keeps the soliton from collapsing. As for the Lee-Wick NTS (see also ) the total energy is expressible as the sum of (i) the volume energy, (ii) the surface energy and (iii) the energy of the fermions:

$$E=\frac{4\pi }{3}V(\varphi _{in})r_b^3+4\pi sr_b^2+\alpha \frac{N_f^{4/3}}{r_b};\alpha \frac{1}{2}(\frac{3}{2})^{5/3}\pi ^{1/3}O(1)$$
$`(54)`$
For a vanishing $`V(\varphi _{in})`$, the minimum of (ii) and (iii) give the total energy $`EN_f^{8/9}(4\pi s)^{1/3}`$, with an energy per fermion $`(4\pi s)^{1/3}/N_f^{1/9}`$. For $`V(\varphi _{in})>0`$, large $`r_b,N_f`$ solutions have energy $`E4(4\pi V)^{1/4}N_f/3`$, with the energy per fermion roughly independent of the number of fermions. Finally, for negative $`V(\varphi )`$, NTS’s exist only for a small fermion number $`(\overline{N}_f)`$ and for small enough $`|V|`$. In this last case, the existence of NTS’s requires the surface term to dominate over the volume term. This is possible only for small $`r_b,\overline{N}_f`$. The energy per fermion is again $`(4\pi s^{})^{1/3}/\overline{N}_f^{1/9}`$. Energy per fermion for a large fermion number configuration is smaller, the closer the minimum of $`|V|`$ is to zero.

This can form the basis for a prescription to sort out the cosmological constant problem. We first recall that in standard GTR, the cosmological consatnt problem is sought to be solved by first making the cosmological constant dynamical: i.e. a constant of motion. This is done by using (for example) Witten’s  prescription by Coleman, Hawking and Weinberg . One adds to the action a term:

$$I_F=\frac{1}{48}d^4x\sqrt{g}F_{\mu \nu \rho \sigma }F^{\mu \nu \rho \sigma }$$

with $`F_{\mu \nu \rho \sigma }`$ the exterior derivative of a three form gauge field $`A_{\nu \rho \sigma }`$:

$$F_{\mu \nu \rho \sigma }=_{[\mu A_{\nu \rho \sigma }]}$$

using the fact that $`F_{\mu \nu \rho \sigma }`$ is totally antisymmetric, it is proportional to the permutation tensor. The field equation for $`A_{\nu \rho \sigma }`$ implies that the proportionality factor is a constant. This reduces the action to:

$$I_F=\frac{1}{2}c^2d^4x\sqrt{g}$$

With $`c`$ a constant of integration. This makes the effective cosmological constant a constant of integration. This being the case, then in a quantum theory we expect the wave function of the universe to be a superposition of states with different values of the cosmological constant. Hawking  has proposed that the application of conditional probability that homosapiens would come into being and study physics in a universe then solves the cosmological constant. However, such an anthropic principle is not a dynamical principle in the first place. Even if the evolution of an intelligent species were related to a dynamic prescription, it would be difficult to implement it in standard GTR where every value of an effective cosmological constant corresponds to a distinct FRW solution that can not be dynamically compared with a solution with a different cosmological constant.

In the structure that we have described, the situation is simpler. One can start with the action eqn(2) and consider NTS’s with a large average fermion density in the universe. Every NTS is a ball that has an effective gravitational constant vanishing in its exterior and having a unique value in its interior determined by the NMC at $`\varphi _{in}`$ where the effective potential is minimized. The effective cosmological constant vanishes in the exterior and the value of $`V(\varphi _{in})`$ contributes to the interior volume energy of the NTS. One can now play the same game as played in standard GTR and introduce $`I_F`$ as before to make $`V(\varphi _{in})`$ dynamical. Consider the wave function of the universe containing a large (conserved) number of fermions. We can expect the state vector of the universe to be a superposition of states with different $`V(\varphi _{in}`$). From the previous considerations, the minimum energy per fermion would be for the state that has $`V(\varphi _{in})`$ vanishing. The minimization of energy is a dynamic principle that may or may not have anything to do with the existence of apes !! With all NTS’s having having common asymptotics for an arbitrary $`V(\varphi _{in})`$, a dynamical solution that uses an energy minimization prescription makes sense and is easy to implement.

### 2(g) Microscopic domains:

It is interesting to note that were we to ensure bounded NMC inside domains as large as hadrons, one would again get an effective gravity model with no cosmological constant problem. Consider for example a theory described by eqn(42) in which the confined fermionic degrees of freedom were associated with baryons. With the NMC diverging at say $`\varphi =0`$ and the fermions confined to hadronic size microscopic $`\varphi 0`$ domains, there would be effective gravitational effects for baryons. A “cosmological constant” as large as $`\mathrm{\Lambda }_{QCD}(.2.3Gev)^4`$ would merely make a small contribution to the mass of a baryon and could be absorbed in its net mass. The vanishing of the covariant divergence of the matter stress energy tensor implies the geodesic law of motion for material particles. The geodesic law in turn implies the equality of the inertial mass and the passive gravitational mass. The gravitating (active) gravitational mass, however, would depend upon the number and binding energy of non-baryonic matter of the source. Considering the relatively large error bars in the determination of the gravitational constant in a typical Cavendish experiment, the change in the active gravitational mass would not be discernible. However, it is claimed  that the equality of the active and passive gravitational masses is established in lunar ranging experiments to a very high accuracy. For this reason (not withstanding our doubts on the claimed accuracy and interpretation of the above experiments) we feel that the large NTS’s described earlier offer a better promise for a problem free gravity theory.

## 3. Features of a toy cosmology:

In the model described above, a divergent NMC over most of the universe provides for a vanishing cosmological constant. Non-trivial solutions for the scalar field provides for gravitating pockets inside compact domains - with the divergent NMC being restricted to their exterior. A value $`N_f10^{75}`$ for $`s(Mev)^3`$ would give a NTS of a size of tens of kilo parsecs. Such an $`N_f`$ is of the same order as the relic background neutrinos / photons in the universe. Thus a fermion species that decouples very early in the universe when that species is still relativistic, would be sufficient to provide $`N_f`$ for gravitating domains as large as a Halo of a typical large structure (galaxy / local group etc.).

With the effective gravitational constant identically vanishing outside a NTS, one could conceive of a cosmological model that starts with a hot big - bang and evolves as a Milne model  almost from its birth. We assume the existence of conditions conducive for the production of above NTS’s with fermions confined to compact domains where the scalar field is non-vanishing. These NTS’s would then evolve by mutually colliding and coalescing. As the NTS’s coalesce, the fermion number adds up, the NTS becoming larger and could evolve to the current distribution of gravitating domains.

On large scales, the universe evolves with the Freidmann - Robertson - Walker scale factor increasing linearly with time: $`a(t)t`$. This has characteristic features: (i) With $`_0^t𝑑t/a(t)`$ unbounded for any $`t>0`$, there is no horizon problem in the theory. (ii) With the expansion parameter not determining any “critical - density” in the model, there is no flatness (fine tuning) problem. (ii) The model is concordant with standard classical cosmological tests, viz.: the number count, angular diameter and the luminosity distance variation with redshift. The first two tests are quite sensitive to models of galactic evolution and for this reason have (of late) fallen into disfavour as reliable indicators of a viable model. However the magnitude - redshift measurements on SN 1A have a great degree of concordance with $`\mathrm{\Omega }_\mathrm{\Lambda }=\mathrm{\Omega }_M=0`$ (). (iii) With the scale factor evolving linearly with time, the Hubble parameter is precisely the inverse of the age t. Thus the age of the universe inferred from a measurement of the Hubble parameter is 1.5 times the age inferred by the same measurement in standard matter dominated model. Such a cosmology promises consistency with an older universe. (iv) The deceleration parameter is predicted to vanish. (v) early universe nucleosynthesis is not ruled out in a Milne universe. It has been shown  that for a baryon entropy ratio $`\eta 5\times 10^9`$, one gets $`24\%`$ of helium-4 and metallicity quite close to that observed in type II stars and low metallicity clouds. The cosmology thus comes with its characteristic predictions and may well be distinguishable by the next generation of experiments.

## 4. Conclusion:

What we have profiled is a plausible program with an exploratory spirit. We had been looking for a framework that could provide for a spatial variation of the effective gravitational constant as a solution to the cosmological constant problem. It was this search that attracted our attention to remarkably encouraging results in the theory of non - topological solitons \[NTS’s\] in the Lee-Wick model. Given a (large), conserved fermion number and a given surface tension of such a soliton, the smaller the magnitude of the interior volume energy density (given by the value of the minimum of the effective potential), the smaller is the energy per fermion of the NTS. If one adds non-minimal coupling, that diverges at a point that we take (without any loss of generality) as $`\varphi =0`$, and makes the interior volume energy dynamical, the minimizimg energy can provide a dynamic prescription to sort out the cosmological constant problem. The effective gravitational constant also approaches a universal value determined by the NMC at the minimum of the effective potential.

We have demonstrated that in a whole class of scalar tensor theories in which the non-minimal coupling diverges and for which the classical effective potential vanishes at some point, classical scalar field condensates can occur as NTS’s. The effective gravitational constant inside all large domains would approach a universal value and the effective cosmological constant would drift to zero. The dynamical tuning of the effective cosmological constant to a small value and the effective gravitational constant to a universal value are compelling features - enough to explore the possibility of raising the toy model described here to the status of a viable cosmology.

Acknowledgements:

This work has benefitted greatly from discussions with Profs. T. W. Kibble, M. Rocek, P.D. Mannheim, L. Parker and G. W. Gibbons.

Figure caption:

Figure 1:The profile of the potential, V and the effective potential, W. The two can be chosen to match for all $`\varphi >0`$ except for $`\varphi `$ near zero.