A Modiﬁed Dynamical Model of Cosmology I Theory

: Wheeler (1964) had formulated Mach’s principle as the boundary condition for general relativistic ﬁeld equations. Here, we use this idea and develop a modiﬁed dynamical model of cosmology based on imposing Neumann boundary condition on cosmological perturbation equations. Then, it is shown that a new term appears in the equation of motion, which leads to a modiﬁed Poisson equation. In addition, a modiﬁed Hubble parameter is derived due to the presence of the new term. Moreover, it is proved that, without a cosmological constant, such a model has a late time-accelerated expansion with an equation of state converging to w < − 1. Also, the luminosity distance in the present model is shown to differ from that of the Λ CDM model at high redshifts. Furthermore, it is found that the adiabatic sound speed squared is positive in radiation-dominated era and then converges to zero at later times. Theoretical implications of the Neumann boundary condition have been discussed, and it is shown that, by ﬁxing the value of the conjugate momentum (under certain conditions), one could derive a similar version of modiﬁed dynamics. In a future work, we will conﬁne the free parameters of the Neumann model based on hype Ia Supernovae, Hubble parameter data, and the age of the oldest stars.


Introduction
The question of the existence of dark matter and dark energy and their implications at different cosmic scales has been the dominant subject of the contemporary cosmology. The main argument in favor of these two components is that their presence would solve some otherwise baffling discrepancies between observations and predictions of general relativity (GR). See References [1,2]. These predictions are mostly based on the zero-and first-order perturbation equations of GR field equations. See References [3][4][5][6][7] for reviews on perturbation method in general relativity.
For example, the presence of the dark matter particles could explain most notably, among other phenomena, the growth of structures in the universe, the stability of gravitating systems, the rotation curves of spiral galaxies, and the missing mass problem in gravitational lensing. See References [1,8,9] for excellent reviews on the essence of these problems. However, it should be noted that, despite tremendous efforts to detect the most popular candidates of dark matter, such as sterile neutrinos, WIMPs (weakly interacting massive particles), and axions, a conclusive evidence for the detection has not yet been reported [10].
On the other hand, the existence of the dark energy is presumed to explain the large-scale structure and evolution of the cosmos [2,11]. The main concerns here are accelerated expansion of the universe (proved by References [12,13] based on supernova observations), baryon acoustic oscillations, temperature anisotropies in the cosmic microwave background (CMB), and the age of the universe. See Reference [11] for a thorough review on theoretical implications and observational evidence of dark energy. Also, see Reference [14], which reviews dark energy models from different perspectives.
To solve GR field equations, the initial/boundary conditions are usually considered to be trivial. See Referece [1] pages 361-363 for a discussion on this issue. The question of initial/boundary condition in GR is closely related to the problem of the boundary term in GR action. At least from a mathematical point of view, it is well understood that, to derive the Einstein equations, one has to add a surface term to the Einstein-Hilbert action [15]. This ensures a well-posed action, meaning that the imposed boundary condition on the action is compatible with the derived field equations. The main issue in adding a surface term to the Einstein-Hilbert action is the quantity which should be fixed on the boundary. On the other hand, Wheeler's [16] interpretation of Mach's principle suggests that, in addition to the energy momentum tensor, one needs also to specify the boundary condition to derive the structure of spacetime. See also Reference [17] for a review on this idea. This view toward solving Einstein field equations is mostly based on the necessity of mathematical and physical well posedness of the theory rather than philosophical grounds of Mach principle, which are well discussed by Reference [18]. In fact, this interpretation of Mach's principle is closely related to the Cauchy problem of GR [19,20]. There are various formulations of Mach's principle in the literature. For example, Bondi and Samuel [21] distinguish eleven statements, Mach0 to Mach10, of Mach's principle; however, Wheeler's formulation is not included. We will discuss some of these statements below. Also, References [17,22] propose similar ideas to Wheeler's in their attempt to find generally covariant integral formulation of GR.
In their work on the foundations of gravitational theories, Thorne, Lee, and Lightman [23] define boundary conditions as those laws which involve only confined variables. As these authors have illustrated, "if one augments the theory by cosmological demands that φ and ψ [the scalar potentials] go to zero at spatial infinity, those demands are boundary conditions". We will discuss this particular case in the following.
Including boundary dynamics while solving GR field equations could lead to interesting results even at the quantum level [24,25]. For example, Park [26] considers the issue of "changing black hole mass through Hawking radiation" and finds that, by imposing Neumann boundary condition, one could solve the problem. See also a recent work by Witten [27], who shows that the usual Dirichlet boundary condition is indeed not elliptic. Thus, in general, this boundary condition does not provide a well-defined perturbation theory. On the other hand, Witten shows that conformal boundary conditions always lead to a sensible perturbation expansion. Therefore, when we solve GR field equations, it is important to check different possibilities of boundary conditions both for the sake of mathematical well posedness of the model and its compatibility with physical expectations.
In this work, we will impose Neumann boundary condition (B.C.) on GR perturbation equations. Also, the modifications in trajectories of massive and massless particles are discussed. The new term in the equation of motion of this model is found to be proportional to cH 0 = 6.59 × 10 −10 m/s 2 in which c is the speed of light and H 0 is the Hubble constant [28,29]. This parameter is also known as the de Sitter scale of acceleration [30,31]. Thus, the model is essentially a modified dynamical (MOD) model which, in many features at galactic scales, shows some similarities with Milgrom's proposal [32][33][34] known as MOND (modified Newtonian dynamics). See Reference [35] for a review on MOND. Some particular aspects of galactic dynamics, such as observed regularities in the properties of dwarf galaxies [36,37], are illustrated within the context of MOND quite better than the cold dark matter ( CDM ) paradigm. The main implication of the appearance of such an acceleration scale in the equation of motion is that, probably, the physics of the local universe might be affected by the global expansion of the cosmos [8]. Many works try to explain the phenomenological success of MOND by deriving particular versions of it from action principle. For example, a recent article by Vagnozzi [38] demonstrates that MOND can be recovered in the low-energy limit of mimetic gravity. The acceleration scale in this approach changes with length scale in such a way that the model is able to explain dark matter both on galactic and cluster scales.
Using our new modified Poisson equation and the averaging procedure of backreaction cosmology [39][40][41][42], in Section 3, a governing modified Hubble parameter of the present model is derived. Then, in Section 4, the main cosmic implications of the new model is surveyed. Among others, we have derived the evolution of the scale factor, Hubble parameter, deceleration parameter, equation of state, and adiabatic sound speed. In a future work, we will present a data analysis of the model based on Supernovae Type Ia (SNIa) and Hubble parameter data (Shenavar and Javidan, in preparation). The basic parameters of the theory could be derived from such analysis.
Some related explanations and elaborations are postponed to the Discussions section and in the appendices. In Section 5.1, we search for the implications of Neumann B.C. from an action-principle point of view. Prior works by Maldacena [43] and by Anastasiou and Olea [44], who propose the equivalence of Einstein and conformal gravity under Neumann B.C., are discussed. Also, we consider the works of References [15,45] that prove the existence of a well-posed Neumann boundary condition. The surface term introduced by References [15,45] has the peculiarity that it vanishes in a four dimensional spacetime. In particular, we have shown that, by assuming some conditions on scalar potentials, our choice of Neumann B.C. in Section 2 could be derived based on the procedure that References [15,45] define. Discussions regarding the averaging problem and Mach's priciple have been presented in the rest of Section 5. In addition, we have provided detailed analysis of the role of geometrodynamic clocks and EPS ( Ehlers, Pirani and Schild [46]) theorem in Appendix B. Moreover, a discussion on the dimensional analysis of the model is presented in Appendix C. Furthermore, we report a similar semi-Newtonian model in Appendix D, which is built based on some well-motivated theoretical and observational assumptions.
In Reference [28], the existence of a cosmological constant Λ was assumed to ensure the existence of a late time epoch of accelerated expansion. Here, however, the Λ term is abandoned and we will show that the modified Hubble parameter derived from imposing a Neumann B.C.-which includes a term proportional to H similar to Veneziano ghost of QCD (Quantum Chromodynamics) [47,48]-would naturally lead to an accelerated expansion of the universe with an equation of state parameter approaching to w < −1 at later cosmic times. The Neumann model proposes a unification of dark matter and dark energy. Such models have been proposed before. For example, one could see the generalized Chaplygin gas [49], k-essence [50], and Bose-Einstein condensation model [51]. We will compare the results of the Neumann model with these models and MOG (modified gravity of Moffat [52]), which is a model of dark matter (plus Λ as dark energy), in the following.

A New Model of Modified Dynamics
In this section, we use the variational method of Weinberg [7], chapter 5, to find the Taylor expansion of GR field equations. One may also see References [3][4][5][6] for precise descriptions of the general perturbation theory of Einstein field equations. To do so, we will first fix the gauge condition to eliminate the spurious gauge degrees of freedom. Then, we will impose the boundary condition on the remaining physical fields.

Fixing the Gauge and the Boundary Conditions
We assume a flat Friedmann-Lemaitre-Robertson-Walker (FLRW) background metric disturbed by the perturbation of the local universe h µν . Under spatial rotation, the symmetric perturbation metric could be decomposed into a scalar h 00 , a three-vector h 0i , and a symmetric two-index spacial tensor h ij . The tensor h ij could itself be decomposed to a trace part and a trace-free part. The resulting field equations, as Reference [7] has reported in chapter 5, would be very complicated. Also, they would contain some unphysical scalar and tensor modes which are correlated with the gauge freedom of the field equations.
This gauge problem could in principle be solved by working with gauge invariant quantities; however, in this way, the issue of imposing a new boundary condition on field equations would not be physically evident. Another method to eliminate the gauge degrees of freedom is by fixing the gauge which we use here. This process could be significantly simplified if the scalar, vector, and tensor modes of the perturbation metric are treated separately. Consider a general coordinate transformation (1) in which ξ µ is of the same order of the perturbation. To distinguish the behaviour of different modes under this transformation, we decompose the spatial part of ξ µ as with the condition that the ξ V i is a divergenceless vector, i.e., ∂ i ξ V i = 0. As Weinberg [7], pages 235-238, shows that the tensor quantities in perturbation metric and perturbation energy-momentum tensor are gauge invariant, there is no need to gauge-fix these modes. On the other hand, the vector ξ V i could be chosen so that some vector degrees of freedom vanish. Lastly, to gauge-fix the scalar modes, there are various possibilities; however, here, we use the Newtonian gauge because, in this gauge, it is more convenient to adopt a new boundary condition. Of course, conversion to other gauges could be done as usual. In Newtonian gauge, two scalar degrees of freedom could be eliminated by choosing ξ S and ξ 0 . Thus, the remaining scalars lead to the next metric (Weinberg [7] page 239) in which R(t) is the scale factor and the potentials Φ and Ψ are functions of spacetime, i.e., Φ = Φ(t, x, y, z) and Ψ = Ψ(t, x, y, z). Our notation is very close to Weinberg's [7], however, as a deviation from Weinberg's notation, we keep a(t) for the dimensionless scale factor a(t) ≡ R(t)/R 0 in the following. One should note that the two scalar potentials are, in principle, distinct degrees of freedom. In fact, from a physical point of view, the potential Φ is a generalization of the Newtonian potential because it specifies the particle acceleration while the potential Ψ is the 3-curvature perturbation of a constant-time surface. While many works consider the two potentials to be the same, the key assumption of this work is that these potentials could be essentially different.
The resulting field equations, as we will see below, would contain second-order spatial derivatives which certainly need boundary conditions to provide a unique solution. Indeed, by imposing a suitable boundary condition, as Wheeler's interpretation of Mach's principle demands, one could uniquely determine the solution to the field equations. The aim of the present work is to provide a new solution for the perturbation equations based on imposing Neumann boundary condition.
Using the metric in Equation (3), it is now straightforward to derive the local and background terms of Christoffel symbols, Riemann, Ricci, and Einstein tensors. These quantities have already been reported by many authors for the full perturbation metric (which includes vector and tensor terms in addition to scalar potentials) [3,4,6,7]. For the sake of brevity, here, we only report Christoffel terms and elements of Einstein tensor which will be needed in following.
The Christoffel symbols, up to first-order terms, are while the components of Einstein tensor are as follows All terms higher than the first order have been dropped in Equations (4) and (5).
It is possible now to find the Einstein field equation, i.e., G µ ν = κT µ ν (10) in which κ = 8πG is a constant, by using Equations (5) and (6). For example, by using i j components of the Einstein equation, one could find which is reported in many classical literature on cosmic perturbation theory. See, for instance, Bardeen [3] and Mukhanov et al. [4]. Equation (11) manifests that the difference between tidal forces due to Φ and Ψ is due to the anisotropic stress π S . In particular, if π S = 0, then the tidal forces However, the field equations are at hand at this point, the appropriate boundary condition is yet to be imposed [16]. See Figure 1, which shows a typical gravitational source (here a galaxy) in a smooth background and its related boundary. The size of the boundary is assumed to be determined by the particle horizon because we assume that the boundary in Figure 1 includes all the objects which have already been in causal contact with the central galaxy which is under consideration [53]. The matter distribution on large scale is assumed to be homogeneous and isotropic; thus, the surface of the boundary, which lies at very large distance, must have the same condition too. The boundary condition will be fixed on this surface in what follows.
The most general solution to Equation (11) could be written as . Similar to other partial differential equations, here too, one must impose an appropriate boundary condition to find the unique solution to Equation (11). In particular, the fact that the presumed boundary in Figure 1 is assumed to be homogeneous and isotropic everywhere sets the first three functions to zero, i.e., A 1 (t, x) = A 2 (t, y) = A 3 (t, z) = 0. This result could also be obtained by arguing that the perturbations must be statistically homogeneous and isotropic. On the other hand, assuming π S = 0, there is no a priori reason for putting c 1 (t) = 0; thus, the gravitational potential Φ and the 3-curvature perturbation Ψ could in principle be different. In fact, it could be shown that a pure time-dependent function c 1 (t) does not change the homogeneity and isotropy of the system under consideration (see Section 5.3 for the details).
In general, because the dynamics of the universe at its largest scale is time dependent, the difference between Φ and Ψ on the boundary could be, in principle, a time-dependent function. However, as we have shown in Appendix A, a time-dependent function c 1 (t) would result in three independent cosmic equations, i.e., two modified Friedmann equations plus a modified conservation equation, such a system of equations is mathematically possible; however, its dynamical behaviour could be completely different in early cosmic times from that of the standard model.
On the other hand, by assuming c 1 as a time-independent factor, one obtains a model for which only two cosmological equations (out of three) are independent. In this way, our analysis of cosmic evolution becomes more straightforward. See Appendix A for details. Also, as we will show in the following, the existence of this term does not alter the "form" of the governing Einstein equations because its space and time derivatives are zero. The appearance of the new solution, however, is due to imposing a new boundary condition. In the same way, the solutions to the geodesic equation will be modified now; however, the form of the geodesic equation is still the same as before.
It could be shown that, after averaging the full metric, the result contains a term proportional to c 1 . However, by a simple redefinition of the scales, one could prove that the averaged metric is the same as Minkowski metric. See Section 5.3 for a detailed discussion.
If one assumes c 1 = 0, then the usual formalism of the standard cosmology could be derived readily. This choice could be interpreted as imposing the Dirichlet boundary condition in solving Equation (11). From our perspective, by imposing the Dirichlet boundary condition on local perturbation equations, one presumes that distant objects enforce no observable effect on local physics, i.e., no detectable result from Mach's principle. On the other hand, a nonzero c 1 could be seen as imposing a Neumann boundary condition [28]. In fact, a nonvanishing potential at the boundary, i.e., c 1 = 0 in Figure 1, represents the effect of the distant stars on local dynamics. Therefore, in the present model, the Mach's principle provides observable effects for the motion of local objects; however, to be precise, this new model is achieved by following Wheeler's work [16] that emphasized the role of the boundary condition in solving Einstein's field equation. See also a discussion in Section 5.2. Figure 1. A particle (dashed curve) is falling into the potential well of a massive object (a galaxy) in a homogeneous and isotropic background. The sphere includes all the particles which have already been in causal contact with the central mass under consideration; thus, the radius of the spherical boundary is indeed the size of the particle horizon. The presence of distant objects produces a nontrivial boundary condition Φ − Ψ = c 1 on the boundary surface, which we will consider here.
From physical point of view, the Neumann B.C. is quite sensible because it implies that, instead of the difference between the potentials (Φ − Ψ), the difference between the tidal forces, i.e., ∂ i ∂ j (Φ − Ψ), vanishes on the boundary.
Chameleon cosmology too, uses a nonvanishing potential as its boundary condition [54]; however, the procedure is quite different in that case. Also, many recent works presume nonzero anisotropic stress, i.e., η = Ψ/Φ = 1, to search for hints of deviations from the ΛCDM model. See, for example, Reference [55] and references therein. However, it should be clear now that the new effect in the present model does not occur due to a nonvanishing anisotropic stress. The appearance of c 1 is because of the presence of matter at the boundary.
Choosing Neumann B. C.
one could rewrite Einstein's field equations as while the components of energy momentum conservation could be found as follows: (14) in which pure time-dependent terms are distinguished with an underbrace. As one could see from Equations (13) and (14), as a result of adding c 1 to the solution of Equation (11), there appears no modification in the zeroth-order (pure time-dependent) GR equations.
In the rest of this section, we will consider the modifications due to Equation (12) on the trajectories of massless and massive particles. This is essential in understanding the present model, and it has been partially discussed before [28]; however, here, we will use the modified Poisson equation introduced in References [29,56] to build a full cosmological description through finding a modified dynamical term of cosmic energy. After finding this term which we will name ρ c 1 , we will derive the modified Friedmann equation and survey its implications in cosmology.

Trajectory of Particles
To find the trajectory of a massless particle, we will use the method of perturbative geodesic expansion introduced by Reference [57]. See also Reference [58] for an extension of the method. To derive the trajectory of a photon, the tangent vector, i.e., v µ = dx µ /dλ, must be null v µ v µ = 0 and geodesic at any order. Assuming k µ and l µ to be the tangent vectors at zeroth and first orders, respectively, the null condition would provide us with while the time and spatial components of the geodesic equation of the photon (at first order) would become in which we have used the Christoffel symbols of Equation (4). Here, the zero-order terms of the geodesic equation have been neglected. Also, spatial vectors have been shown with an arrow k. The evolution of l 0 , i.e., Equation (17), essentially illustrates the change in the redshift of photons when they pass through a gravitational well. This equation is the basis for the discussion of the Sachs-Wolfe effect in CMB physics. The changes that the presence of a Neumann constant c 1 introduces to the CMB analysis will be discussed in more details in a future work.
The solution to Equation (18), on the other hand, describes the deviation in photon trajectory while passing through a gravitational potential which results in the lensing equation. In fact, using Equations (15) and (16) to replace k 0 and l 0 , respectively, in Equation (18) and neglecting terms like HΦ and H l because of their minute role at local physics, one can finally derive the governing equation of the spatial components as (19) in which the gradient transverse to the trajectory ∇ ⊥ is defined as follows The deflection angle α ≡ − ∆ l k 0 , in which ∆ l = d l dλ dλ, could be simply derived as reported in Reference [28].
We stress the fact that the above modified trajectory is derived based on the method of Reference [57]. However, if one hastily applies the usual lens formula of α ∝ ∇ ⊥ (Φ + Ψ)dλ and replaces Ψ from Equation (12), then there remains no modified term proportional to c 1 . On the other hand, by first deriving the deflection value d l/dλ at any point of the trajectory as Equation (19) and by then integrating over photon geodesic d l dλ dλ, one could retain the role of the boundary condition c 1 , i.e., the effect of the distant stars is preserved. The difference between the two methods is that, based on the method of Reference [57], which we used here, one has to use Equation (16) to replace the value of l 0 into Equation (18), which is a more appropriate approach.
The appearance of the acceleration term c 1 cṘ R in Equation (19) is the result of the coupling between the background and the perturbed metrics. Here, we have recovered the velocity of light c to estimate the order of magnitude of this term. From physical point of view, this term indicates that the evolution of very far objects affects the local physics of gravitating systems. The magnitude of this acceleration is of the order of 10 −10 m/s 2 to 10 −11 m/s 2 depending on the value of the Neumann parameter c 1 . A similar quantity has been repeatedly observed in MOND phenomenology [32,35].
The aforementioned modified lens equation has been tested against a sample of ten strong lensing systems, and it has been shown that the estimated masses of these systems are within the observed values except for the case of the Q0142-100 lensing system [28]. The evaluated mass of this system shows about 7.5% deviation compared to the lower observed bound. This case is not in strong contradiction with the model regarding the possible sources of uncertainties, such as error in the position of the image, the impact parameter, and approximating the lenses as spherical systems.
The above method also changes the geodesics of massive particles. The key point here is that, in GR, one has to use geometrodynamic clocks to measure spacetime intervals and physical observables such as velocity, acceleration etc. [28]. Geometrodynamic clocks use particles in free fall and light signals to measure space and time. See Reference [59], chapter 5, for a review on spacetime measurement through geometrodynamic clocks. Thus, a modification in the light's geodesics would consequently result in a change in measurement's outcome. See Appendix B for a related discussion on geometrodynamic clocks and the implications of the key EPS theorem [46]. Also, it is possible to generalize the equation of motion of a single particle to a system of particles. See "Section 3-Systems of Particles" of Reference [29] for the details. Then, the appropriate potential of MOD is derived, and it has been shown that this model is governed by the next modified Poisson equation: where plays the role of the missing mass at galactic scales. Appendix B provides a thorough derivation of Equation (20). The parameter c in Equations (20) and (21) is recovered to manifestly show that δρ c 1 and δρ m have the same dimension. At galactic scales, the time evolution of H(t) and δρ m (t, x ) is usually considered negligible because the time period that one considers a typical galaxy is much shorter than the cosmic time for which H(t), and δρ m (t, x ) changes significantly. Of course, the problem of galactic evolution is an exception. Also, in this work, we deal with cosmic evolution; thus, it is important to include the time evolution of these factors. We note that δρ c1 is in fact determined by δρ m and that it is actually a new energy between particles; therefore, δρ c1 has no independent existence apart from the matter distribution. Strictly speaking, δρ c1 does not clump; only matter density δρ m clumps. Therefore, no interaction between δρ c1 and δρ m would exist. However, we will continue to show this energy as δρ c1 because, in this way, it is more convenient to incorporate it within the field equations and to derive the modified Friedmann equations.
One could evaluate δρ c 1 by integrating Equation (21) over the space in which the matter exists. This is a nontrivial work; however, it could be done at least numerically. The important point is that the integral in Equation (21) has the dimension of mass per scale length ∝ M/R p in which R p is the particle horizon and M is the mass within it. See Figure 1. Thus, the mass M could be simplified in Equation (21). On the other hand, if we show the current size of the particle horizon as L 0 and its time evolution as the dimensionless function l(t), i.e., R p = L 0 l(t), then the modified density δρ c 1 could be written as which represents the effect of distant cosmic objects on local physics. We assume for the sake of simplicity that the constant dimensionless parameter which emerges from the integration in Equation (21) has been absorbed in the parameter L 0 . The logic behind this is the fact that the parameter L 0 will disappear when we write the equations based on dimensionless densities. If the time evolution of the scale factor was known at this stage, then we could determine the size of the particle horizon from R p = L 0 l(t) ≡ c da/aȧ. However, we still have some clues about the past epochs of the cosmic evolution. For example, when the universe is in the phase of radiation domination (or matter domination), the share of the new term ought to be negligible. Thus, as we will show in Section 4, the evolution of the scale factor could be found as a(t) ∝ t 1/2 and a(t) ∝ t 2/3 for radiation and matter-dominated universes, respectively. Therefore, the particle horizon could be calculated as R p = 2ct for the former and R p = 3ct for the latter.
Although we could not determine the parameter l(t) in recent cosmic era, the important point regarding l(t) is that this is a function of time. For an expanding universe which starts with a singularity and has no contraction, it could be easily understood that the scale factor is a reversible function of the cosmic time. Therefore, the parameter l(t) could be written as a function of the dimensionless scale factor a instead of t. In this work, we model the ratio ζ ≡ 1/l(t) in δρ c 1 based on the next power-law function of the scale factor ζ = a ε . (22) in which the value of ε would be determined through data analysis. This toy model proves to be useful especially for low-redshift observations. A preliminary data analysis of the model, with the above assumption, has been reported in the first Arxiv version of this work. See "Section 5-Data Analysis", especially Tables 1 and 2 and Figures 9-14, of arXiv:1810.05001v1 for the detailed discussion. However, we should mention that, for a problem with a broader range of redshift, for example CMB physics, one might need a more general function of the scale factor to describe ζ.
Using the above definition of ζ, the modified energy density of this model is as follows: In a future work and based on the Hubble parameter and SNIa data, we will report that the value of ε is most probably about 0.15 (Shenavar and Javidan, in preparation). The effect of changing ε on data fittings will be also discussed. On the other hand, as mentioned above, the exact value of L 0 will not be needed in building a cosmological model because the dimensionless density parameter related to δρ c 1 is independent of the scale length L 0 . In classical hydrodynamics of fluids, within the context of Newtonian mechanics, one does not expect that adding a constant to the gravitational potential result in any change in the final physics. The reason is that the gravitational effects are introduced to the governing equations through the force (not the potential). On the other hand, in general relativity, the two scalar potentials are basically independent. Regarding this point, the coupling between zero-order and first-order terms of metric, i.e., H and c 1 , respectively, introduces a nontrivial term of the order of de Sitter scale of acceleration cH 0 as we saw above. Finally, it is easy to observe that the MOD model reduces to old results in the limit of Dirichlet B.C. c 1 = 0; however, the vice versa is neither correct nor necessary according to "correspondence principle". Although adding a constant to the potential in Newtonian gravity would provide no observable effect, it is still possible to build a semi-Newtonian model of modified dynamics based on some simple and well-motivated assumptions. See Appendix D for the details. The result would be almost the same; however, the main difference is that, in this latter approach, the value of the new acceleration should be determined from observations.

Modified Hubble Parameter
In this section, we want to find the average of the perturbations to derive the governing equations at cosmic scales, i.e., the Friedmann equations. Literature related to backreaction cosmology have produced a rigorous mathematical framework to study different perturbative orders and their averaging process. This is expected because backreaction cosmology tries to find the effects of inhomogeneities on background metric to explain the accelerated expansion of the universe [39][40][41][42]. Thus, here, we will use the averaging method usually presumed in backreaction cosmology. In the following, the local and background quantities are displayed by l and b indices, respectively. Assuming this, the local perturbation equations can be written as (24) which, by considering the negligibility of background quantities such as ρ, P, H 2 ,R, compared to the local density of galaxies and clusters, could be rewritten as in local universe. See Table 5 of Reference [60] for order of magnitude estimations of the ratios of local-to-global density and pressure at different scales.
On the other hand, the time evolution of the cosmos, at its largest scales, could be derived by finding the spatial average of the inhomogeneities. However, the process of spatial averaging of Einstein field equations is still debatable [61,62], in the literature related to cosmology (especially backreaction cosmology), it is customary to use the next definition to find the spatial mean of any inhomogeneities A(t, x) [63]: where h is the determinant of the perturbed metric on hypersurface of constant time. In this way, the background Einstein equation is as follows by which one could see that because the spatial average of all space-dependent components of the Einstein field equations-including Φ, δρ m , δP, δu, and π S -would vanish; only pure time-dependent terms would survive. These terms are shown with underbrace in Equations (13) and (14). Also, the term δρ c 1 would survive because, as one could see from Equation (23), this term is pure time dependent. We will show the average of δρ c 1 as ρ c 1 , i.e., ρ c 1 ≡ δρ c 1 . Using this process of spatial averaging, one could derive the following from Equation (13) and (14): in which now the total energy of the universe consists of radiation ρ r , baryonic matter ρ m , and the term due to our modification ρ c 1 which is defined in Equation (23). In other words, we have In the present work, we neglect the cosmological constant Λ in contrast to the procedure that was used in Reference [28]. The reason is that, as we will discuss below, a term proportional to H in the Friedmann equation could provide a deceleration parameter about −1 in recent cosmic time. See Section 4.3 below.
The first two equations in Equation (27) can be rewritten as in which we have resumed the role of the curvature constant k for the sake of completeness and future references. Here, the parameter k is different from the zero-order wave vector k µ (and k) in the previous section. The values of k = 1, 0, −1 refer to closed, flat, and open universes, respectively. These equations have essentially the form of Friedmann equations because, as we saw above, the form of Einstein equations remains unchanged in our approach. However, the energy content of the universe is now changed due to the presence of ρ c 1 , i.e., the modified dynamical energy. In addition, one can rearrange the conservation equation, i.e., the third equation of Equation (27), as which could be dealt with more easily because the time dependency of the parameters is now implicit.
In this work, we assume that the energy density of the universe is only due to radiation ρ r , baryonic matter ρ m , and Neuamnn term ρ c 1 . We will find that, in this modified cosmological model, a large portion of the energy content of the universes is due to this last term. By including the effect of Neumann B.C. through ρ c 1 , one has the total density as ρ = ρ m + ρ r + ρ c 1 . The cosmological energy components are considered to be noninteracting; thus, the conservation equations, i.e., Equation (29), holds separately for each component. Also, one could assume an equation of state of the form P i = w i ρ i in which we have used indices i = m, r, c 1 for matter, radiation, and c 1 components respectively. Then, it is possible to show that, for dust with P = 0, one has ρ ∝ R −3 while, for radiation with P = 1/3ρ, one derives ρ ∝ R −4 . The total equation of state parameter, i.e., w, is time dependent. Using the above results and the expression of ρ c 1 = c 1 c πGL 0 Hζ, one could find the evolution of the total density as in which ρ 0,m and ρ 0,r are respectively the matter and radiation densities at the present time. Knowing that ζ is dimensionless, one could easily check that the last term on the rhs has the dimension of density.
It is worth noticing that there is a fundamental difference between ρ r and ρ m on one side and ρ c 1 on the other. The point is that, although ρ r and ρ m are related to physical substances, namely radiation and matter, the term ρ c 1 is due to the modification of dynamics. In fact, as one can see from the definition of ρ c 1 , this term is related to Hubble parameter H and scale factor; thus, from conceptual point of view, it might be even more appropriate to put it on the left-hand side of the Friedmann Equation (28) because it is more similar to H 2 (In this respect, the present model is within the larger class of dynamical dark energy models.). However, since we want to compare our results with those of ΛCDM, we treat this term as a density and keep it on the left-hand side of Equation (28). Of course, the results of the two approaches would be identical because they are mathematically equivalent.
In cosmology, it is well-known that working with dimensionless densities, i.e., the density parameters Ω, is easier than working with the physical densities ρ i . Defining the density parameters as one could rewrite the first Friedmann equation in Equation (28) as follows: is the normalized scale factor. In definition of a(t), the parameter z represents the redshift. From Equation (31), one could see that, for the curvature density parameter, we have While other density parameters are positive by definition, the curvature density parameter could be positive, negative, or zero depending to the the value of ∑ i Ω i (t). Equation (31) could be simply solved to derive the Hubble parameter H(t). By doing so, one finds which is the key formula for our following discussions. Imagine that Ω 0,k 0, i.e., k 0. Then, the plus sign in Equation (32) represents an expanding universe, i.e., H > 0, while the minus sign shows a contracting cosmos with H < 0. On the other hand, when Ω 0,k < 0, i.e., k > 0, then at some normalized scale factor a, the value under the square root in Equation (32) could be negative, which makes the Hubble parameter imaginary and so the real part of the scale factor would be oscillatory. However, when this value is nonnegative, then the Hubble parameter could be positive or negative, depending to the sign of ± in Equation (32), and so we might have an expanding or contracting universe in this case. As we will see in the next section, Ω k can be neglected in our universe.
It is insightful to compare the behavior of the current model with that of the ΛCDM model. The Hubble parameter in ΛCDM cosmology is as follows in which Ω 0,Λ ≡ Λ/3H 2 0 is the density parameter of the cosmological constant while the definitions of the other parameters have been presented above. Again, it is possible to introduce an expanding (contracting) universe with H ΛCDM > 0 (H ΛCDM < 0). Of course, we will consider an expanding universe.
As stated above, in Neumann cosmology, the term responsible for the accelerated expansion changes with time. Models of the universe with variable cosmological terms, i.e., dynamical dark energy models, have been repeatedly proposed before. See Reference [64], especially Table I, for a review on a variety of these models. In fact, there have been previous papers suggesting models with a dark component of which its energy density is proportional to Hubble parameter, i.e., Ω ∝ H(t).
The behavior of such models are somehow similar to the present model, though their origin is far different. These models are proposed based on Veneziano ghost of QCD. This ghost is considered as unphysical in the usual Minkowski spacetime. However, it could lead to interesting cosmological consequences in dynamical spacetimes or spacetimes with nontrivial topology. See References [47,48] and references therein for more details. Assuming Λ QCD ≈ 100 MeV as the mass scale of QCD, one could show that Λ 3 QCD H is of the order of observed dark energy density. This coincidence makes the Veneziano ghost model particularly interesting. See also a relevant discussion in Appendix C.
Cai et al. [47] have shown that Veneziano ghost leads to a late time de Sitter phase of the universe and its accelerated expansion at z ≈ 0.6. They have used observational data of big bang nucleosynthesis, baryon acoustic oscillation, type Ia supernovea, Hubble parameter data, and cosmic microwave background to estimate the free parameters of the model. The best-fit values of the parameters of Veneziano dark energy model lead to a minimum of the probability distributions χ 2 , which is about 10% larger than the minimum of χ 2 derived from ΛCDM. In both models. the existence of the dark matter is presumed. Their study has been expanded by Reference [48] to include a more general model which includes two terms proportional to H and H 2 . This latter model too approaches a de Sitter phase while it begins to accelerate at z ≈ 0.75. Also, according to Reference [48], the data analysis of the growth factor based on this latter model shows a fit as well as those of standard model.
It is worth mentioning that, by imposing Dirichlet B.C., i.e., simply putting Ω 0,c 1 = 0 in our equations, Equations (31) and (32) reproduce the Hubble parameter of the conventional cosmology without a cosmological constant. Therefore, our model includes ΛCDM without Λ. However, as we will see in the following, this model could provide an accelerated expansion without assuming the cosmological constant.

Cosmic Implications of the New Model
In this section, we will outline the main cosmological consequences of the present model. Throughout this section, we will assume that the universe is flat, i.e., Ω 0,k = 0. Also, in all plots of this section, the matter density of the MOD model is assumed to be Ω 0,m = 0. 15 and Ω 0,m = 0.25. These values are the lower and upper limits of matter density according to our data analysis in Shenavar and Javidan (in preparation). Also, we assume ε = 0.15.
For the sake of comparison, the results of ΛCDM model is also provided for which we have assumed the total mass density, the baryonic plus CDM densities, as Ω 0,b+cdm = 0.3. Here, for both models, the present radiation density is considered as Ω 0,r = 5 × 10 −5 .
We first start by discussing the evolution of the scale factor. Then, the behavior of the Hubble parameter is surveyed and it is proved that the deceleration parameter converges to < −1 for a flat universe. This behavior ensures an accelerating expansion at later cosmic times. Moreover, we will derive the equation of state parameter and the sound speed. Furthermore, it is shown that our model reproduces a viable sequence of cosmic eras. Also, the luminosity and angular diameter distances are found to be significantly different at z 2, making a good opportunity to contrast the model with ΛCDM at high redshifts in future works.

Evolution of the Scale Factor
In a flat expanding universe, one could see from Equation (32) that the evolution of the dimensionless scale factor is governed by in which τ ≡ H 0 t is the dimensionless time. This differential equation is solved numerically for two values of Ω 0,m and the initial value of a(0) = 1, and the results are reported in Figure 2. The evolution of the scale factor in the ΛCDM model is also shown by a solid line which closely follows a Neumann cosmology with 0.15 Ω 0,m 0.25. Apart from this numerical solution, one could find the time dependency of the scale factor in some special cases. For example, in early universe, the share of the radiation dominates the right hand side of Equation (34). In this case, the scale factor could be derived as a(τ) ∝ τ 1/2 . Also, in later times, the matter term dominates the rhs of Equation (34) and so we derive a(τ) ∝ τ 2/3 . Finally, at later epochs, the terms corresponding to matter and radiation vanish and Ω c 1 start to govern the cosmic evolution.
Here, if we have ε = 0, then the solution to Equation (34) would be exponential, i.e., a(τ) ∝ exp(Ω 0,c 1 τ). These last three statements could be easily derived from Equation (34). On the other hand, in the case of a nonzero ε, we need to numerically solve Equation (34) to derive the scale factor.

Hubble and Deceleration Parameter
From Equation (32), one can readily plot the Hubble parameter as a function of the scale factor as you may see in Figure 3. In early cosmic time, the behavior of Hubble parameter in MOD is quite similar to that of the ΛCDM model, though the value of Hubble parameter in the former is generally smaller than the latter for a < 1. However, the behavior of H for a > 1 is completely dependent to the future evolution of ζ (or ε). The curves in Figure 3 are graphed for ε = 0.15; however, in the case of ε = 0, we would have a flat Hubble parameter for a > 1. The deceleration parameter q is defined as The Friedmann equations have been used to derive this last expression. The evolution of q is graphed in Figure 4 in which one observes q = 1 in radiationdominated era, then q = 1/2 in matter-dominated epoch, and finally q = −1 − ε in future cosmic times. Also, to plot q as a function of the normalized scale factor, we have used the next identity: which converts time derivative of any quantity to its derivative with respect to a. The negativity of the deceleration parameter at present cosmic time is indeed needed to fit observational data to the model as is discussed by Shenavar and Javidan (in preparation). It is also worth noting that, for a model with ε = 0 (similar to QCD ghost model), one would have a similar behavior except that, in future cosmic time, one derives q = −1. This is not plotted in Figure 4. The negativity of the deceleration parameter in late cosmic time has an important effect on the arguments related to the value of the current cosmic curvature. In fact, using the definition of Ω k , one could derive the time derivative of the curvature density parameter asΩ k = 2Ω k Hq. Then, if Ω k has some nonzero value at a very early universe, the late-time negativity of q ensures that the value of Ω k becomes negligible at present. However, by assuming c 1 = 0, i.e., imposing the Dirichlet boundary condition, the parameter q would always be positive and the Ω k would not be negligible. This is one of the reasons that one needs cosmological constant Λ in the standard model.

Evolution of the Equation of State Parameter and the Sound Speed
In a flat universe, the total equation of state parameter w t could be derived from Friedmann equations as follows: Now one could simply plot w t (a) as shown in Figure 5 by using the conversion formula of Equation (36). As this graph shows, the value of w t is initially very close to 1/3 because the universe starts from a radiation dominated era. Then, in a matter-dominated era, w t approaches very rapidly to zero. Finally, at later times, w t (a) converges (from above) to w = −1 − 2ε/3, which results in an accelerating expansion of the universe. Therefore, in a completely Ω c 1 -dominated universe, the equation of state must be less than −1 for ε > 0. In this regard, the present model is similar to phantom models of dark energy; however, the existence of a big rip in future, which is the characteristic sign of a phantom scenario, crucially depends on the value of ε at that time. In fact, it is more likely that big rip never happens in Neumann cosmology because all forces of this model are attractive. See the modified Poisson Equation (20) above. In contrast, in a phantom model, the big rip happens when the repulsive force of dark energy rips apart all structures down to subatomic particles [65]. See also Reference [14] for a thorough review. The issue of the behavior of ε and, thus, the possibility of a doomsday scenario in future must be resolved with a more careful study to determine the behavior of ζ(t) = 1/l(t). See Vagnozzi et al. [66], who compare dynamical dark energy models with w ≥ −1 and w < −1. The authors show that, unlike models which allow values of w < −1, a dark energy component with w ≥ −1 is unable to reduce the tension between observables of high redshift and direct measurements of H 0 . We will investigate this matter in Shenavar and Javidan (in preparation). It is worthy to note that the curves with smaller Ω 0,m tend to w < −1 more quickly because the share of Ω 0,c 1 is larger in these cases. Also, we point out that the parameter w t of the ΛCDM model coincides (approximately) with 0.15 Ω 0,m 0.25 of the present model. See Figure 5. Another important quantity to be compared with its ΛCDM counterpart is the sound speed. For a barotropic fluid, in which the pressure is only dependent to the density ρ, the sound speed is defined as c 2 s =Ṗ/ρ. This quantity is sometimes called adiabatic sound speed because, in such medium, the entropy per particle is assumed to be constant. As mentioned above, the only physical quantities in MOD model, as in ΛCDM, are radiation and (baryonic) matter. Thus, the sound speed could be found from Friedmann equations as c 2 s = (Ṗ r +Ṗ m )/(ρ r +ρ m ). Using this equation and changing the time derivatives to derivatives with respect to scale factor Equation (36), one could plot the sound speed as Figure 6. The sound speed squared in a radiation-dominated era is about 1/3 while, in a matter-dominated era and beyond, tends to zero. The important point here is the positivity of the sound speed of MOD which results in a bounded (non-exponential) rate of the structure formation. Also, from c 2 s , one could find the angular size of the sound horizon of MOD and compare it with that of the ΛCDM model. This will be done in future.

Sequence of the Cosmological Epochs
It is easy to check that the present model starts from an unstable radiation-dominated mode, then proceeds toward a matter dominated era, which is also unstable, and then ends its evolution at a stable epoch of Ω c 1 domination. See Figure 7. In this plot, we have assumed that Ω 0,r = 5 × 10 −5 , Ω 0,m = 0.25, and Ω 0,c 1 = 1 − Ω 0,m − Ω 0,r . One could easily check that, by choosing other fractions of Ω 0,c 1 and Ω 0,m and by keeping the radiation density at the order of Ω 0,r ≈ 10 −5 , similar behaviors would emerge. In some alternative cosmologies, the sequence of cosmological epochs is not necessarily as described above. If this happens, one might encounter some problems in dealing with the formation of structures in a matter-dominated era. For example, Amendola et al. [67] have proved that f (R) models of the forms f (R) = αR −n and f (R) = R + αR −n do not show a viable matter-dominated epoch prior to a late-time acceleration for any n > 0 and n < −1. In some cases of f (R) theories, the matter dominated epoch is replaced by an era of cosmic expansion in which the scale factor varies as a ∝ t 1/2 , which is cosmologically unacceptable.
Moreover, in the case of scalar-tensor-vector theory of Moffat [52], Jamali and Roshan [68] have shown that there are two radiation-dominated eras, two matter-dominated epochs, and two late time-accelerated phases. In the matter-dominated phases, the growth of the scale factor is found to be as a(t) ∝ t 0.46 and a(t) ∝ t 0.52 slower than the growth in standard model and the present model reported above as a(t) ∝ t 2/3 . However, it should be mentioned that, as [69] has reported, the standard MOG possesses a valid sequence of standard cosmological eras.
Because the dimensionless density parameter of the Neumann term is inversely proportional to Hubble parameter, i.e., Ω c 1 ∝ 1/H(t), it plays no significant role in an early universe which is dominated by radiation and then matter. A similar situation exists in the context of the standard model of cosmology, in which Ω Λ is negligible until recent time. However, in the case of ΛCDM, one deals with the cosmological constant which is very tiny Λ ≈ 10 −52 m −2 , while the present model is built on introducing a dimensionless parameter with a value much closer to unity c 1 = 0.065. Also, the density of the cosmological constant, ρ Λ = Λc 2 /8πG is always a constant while the density of ρ c 1 varies with time according to Equation (23).

Luminosity and Angular Diameter Distances
To compare our model with observational data, we need to measure distances in the universe. Shenavar and Javidan (in preparation) study the behavior of luminosity distance while the angular diameter distance is studied by Reference [70]. Both distances are defined for a flat universe. The plot of luminosity distance D L and angular diameter distance D A for two values of Ω 0,m are shown in Figures 8 and 9, respectively. The analogous curves for ΛCDM model too are presented by a solid line in both figures.  As it is shown in the left panel of Figure 8, the luminosity distance of the present model converges to a constant value at z 10 4 while the same curve for the ΛCDM model varies linearly. Please note that the left panel is logarithmic on both axes. Fortunately however, to compare the present model with the standard one, it is not needed to observe objects at very high redshifts. The reason is that, as it is shown in the right panel of Figure 8, the luminosity distance in the present model begins to differ significantly from that of the standard model at z ≈ 4. If the baryonic content of the universe at the present time is larger than a minimal value, say Ω 0,m 0.15, then the deviation between the two models might be detectable with a precise data analysis at even z ≈ 2.0. Therefore, it seems that the luminosity distance provides one of the best possibilities to check the viability of the present model. In conclusion, we should say that, although the luminosity distance of the present model is rising with increasing z, it is always smaller than D L of ΛCDM at the same redshift. In fact, D L of ΛCDM acts as an upper boundary of D L predicted by the present model.
One of the most important consequences of changing the value of D L is that now the absolute visual magnitude M V of the objects observed at high redshifts has to be modified (compared to the values derived from standard model). Assuming the apparent visual magnitude as m V , the relation between these two magnitudes is in which D L,0 is some constant distance, say 10 pc for near objects or 1 Mpc in cosmological measurements. Then, one could simply prove that in which M V,MOD is the absolute visual magnitude in our modified model, M V,ΛCDM is the absolute visual magnitude in the standard model, D L | MOD is the luminosity distance in the present model, and finally D L | ΛCDM is the luminosity distance of the standard model. As it is clear from above equation, the fact that the luminosity distance in the present model is approximately equal to or smaller than that of the standard model leads to the fact that the absolute visual magnitude of the objects would be approximately equal to or larger than the absolute visual magnitude of the standard model. This point must be included especially when one deals with objects at higher redshifts. In addition, the change in D L could affect the mass estimations in lensing problem if they lie at relatively high redshifts. In the lensing sample which was used in Reference [28], the lenses lie at redshifts 0.2 ≤ z ≤ 0.86, which are relatively low redshifts. Depending to the values of Ω 0,m and Ω 0,c 1 , the mass estimation of lenses at higher redshifts might be affected.
The angular diameter distance of the present model too is (almost) always smaller than D A of ΛCDM, as shown in Figure 9. The value of D A for both models decreases when z 2, though the decline of D A in the present model occurs more rapidly. The angular diameter distance provides a good opportunity to confront different cosmological models through the famous "angular size-redshift" problem. See Reference [71] for ane xample. A comparison between ΛCDM model and some other cosmologies, including the present one, based on angular size-redshift data is provided by Reference [70].

Discussion on the Boundary Condition of GR
Assuming Neumann boundary conditions to solve Equation (11), which results in the solution of Equation (12), is valid as a mathematical possibility. However, this choice of boundary condition would be more motivated if we could enforce it based on some physical backgrounds, i.e., action principle, for example. In fact, the debate on the proper boundary condition of Einstein-Hilbert action dates back to early days of the introduction of general relativity. See Reference [72] for a review on this subject. Especially, one could mention the discussion between de Sitter and Einstein on this matter, which could be considered as the root of relativistic cosmology. In the heart of the debate lies the possibility of a static universe, presumed by Einstein before Hubble's discovery of the expansion of cosmos [1]. To achieve the relativity of inertia, Einstein proposed a metric of the form g 0i ≡ ∞ at spatial infinity while other components are zero. However, de Sitter criticized this point of view because it leads to the notion of Newtonian "absolute time" and some invisible masses. In the context of GR, this could be considered as the first mentioning of "dark matter". In fact, through trying to solve the problem of boundary conditions at spatial infinity, de Sitter and Einstein introduced two distinct cosmic models, the first two models of relativistic cosmology [72].
From a pure mathematical point of view too, the issue of boundary condition seems to be quite controversial in GR. In fact, since Einstein's field equations contain second-order derivatives of the metric tensor, it could be shown that the Einstein-Hilbert action is not well posed, i.e., the boundary condition of this action is not compatible with the obtained field equations [15]. To cure this problem, however, one could eliminate the surface terms by adding a boundary term to Einstein-Hilbert action. For example, York [73] and Gibbons and Hawking [74] proposed the next action with a surface term which is invariant under diffeomorphism. In this action, known as YGH action, R is the Ricci scalar, h is the induced metric on the boundary ∂M with the coordinates y i , K is the extrinsic curvature of the boundary, and g is the determinant of the metric tensor. Also, we have γ = +1 for time-like and γ = −1 for space-like boundaries. This action presumes Dirichlet boundary condition by killing all the normal derivatives of the metric tensor on the surface. See Reference [15,45] for a review. It is also worth mentioning that, according to Reference [75], there could be infinitely many boundary terms to make the above action well posed. Therefore, it is not uniquely determined.
On the other hand, by assuming Neumann B.C., Chakraborty [15] and Krishnan and Raju [45] have shown that the action takes the following form: with the unique property that the surface term vanishes in D = 4 dimension. Krishnan and Raju [45] suggest that this property of S N could be interpreted as follows: "Standard Einstein-Hilbert gravity in four dimensions, without boundary terms, has an interpretation as a Neumann problem". In addition, Equation (42) is interesting from the point of view of gravitational theories in higher dimensions since it singles out the D = 4 dimension. In other words, our four dimensional spacetime is not merely one of the possibilities; i.e., D = 4 is the dimension that one has the elegant Einstein-Hilbert action (without any surface term) if one imposes the Neumann boundary condition. Krishnan and Raju [45] have also shown that, in a three-dimensional space, the Neumann action S N becomes the Chern-Simons action, which is another interesting feature of this action. In addition, Maldacena [43] has proved that, by using Neumann B.C., one can derive the semiclassical (or tree level) wavefunction of the universe in 4-D asymptotically de-Sitter or Euclidean anti-de Sitter spacetimes. Since, conformal gravity has many solutions, Neumann B.C. seems to select the Einstein solution out of all possible choices. A more elaborate derivation of the equivalence between Einstein theory and conformal gravity by assuming Neumann B.C. is also provided by Reference [44].
Chakraborty [15] shows that, to impose the Neumann B.C., one has to fix the momentum conjugate on the boundary: In what follows in this subsection, we will impose the Neumann B.C. on cosmic perturbation equations and then derive the modifications due to this new boundary condition on cosmic equations. To do so, assume a metric of the form Please note that, unlike metric (Equation (3)), the scale factor R(t) is not included in the potential ψ. See Reference [60] for a discussion on the difference between the two metrics. For this metric, the momentum conjugate could be derived as in which the first term on the rhs is of zeroth order while the rest represent the perturbed terms. Neumann boundary condition must be imposed on physical degrees of freedom, i.e., the potentials φ and ψ. Also, presuming that the time derivative of the 3-curvature perturbation ψ is negligible, i.e.,ψ/H ≈ 0, we would have for the local conjugate momentum. Remembering that ψ = R 2 Ψ and assuming l Π ij = 2Ṙc 1 at any time, one would derive Φ − Ψ = c 1 as suggested in Equation (12) (neglecting the anisotropic stress). We note that the Dirichlet B.C. could be reproduced simply by replacing c 1 = 0 in Equation (12). In deriving the above condition from the Neumann boundary of Equation (43), we assumed the negligibility of the time variation of ψ. The validity of this presumption at different scales needs to be studied more carefully and we do not consider it here. Also, another critical assumption is that we only imposed Neumann B.C. on potentials (perturbations) and not the background term, i.e., −2Ṙ, which is related to the geometry of the space. This is justified because we prefer to maintain the underlying geometry intact; otherwise, the whole method of perturbation theory which we used here needs to be revised. In conclusion, the above method provides a possible interpretation for the Neumann boundary condition from an action-principle point of view. The other possible interpretation which was mentioned above and is based on the equivalence between Einstein and conformal gravity would not be discussed here. We should point out that boundary conditions, such as mixed boundary condition, are also possible. See, for instance, Peebles et al. [76] and Peebles [77] impose mixed boundary condition to study the dynamics of the local group.
It is interesting that the possibility of a Neumann B.C. to solve GR field equations has been rarely discussed before. One of the reasons for this shortcoming is that proving the well posedness of GR field equations under a particular boundary condition has been found to be a very difficult task. The well posedness, i.e., the existence and uniqueness of the solution in a small neighborhood, of GR under Cauchy B.C. has been proved by Choquet-Bruhat for the first time and after a long debate. A review on this proof which is based on harmonic coordinates could be found in Reference [19]. Then, Choquet-Bruhat and Geroch [20] proved the theorem of the global existence and uniqueness of GR. See also Reference [78] for a good introduction to this matter. As the authors of Reference [79] have discussed too, the Cauchy problem starts with assuming (h ij , Π ij ) as a complete set of Cauchy data (initial data). It is not yet clear that, if we change the boundary condition, the proof of "the existence and uniqueness of solutions" would remain intact or not. However, a systematic treatment is needed to prove (or disprove) the well posedness of GR with a Neumann B.C., which is beyond the scopes of the present work.
Another reason for the usual discarding of Neumann B.C. is that the meaning of this boundary condition in a four-dimensional spacetime is not necessarily evident. For example, in the field of numerical relativity, there have been some attempts to impose Neumann B.C. to prevent the violation of constraints [80] or to investigate the numerical stability of Cauchy evolution [81] under the new boundary condition. To impose Neumann B.C., Kidder et al. [80] placed restrictions on the normal derivatives of some characteristic fields while [81] determined ∂ z Φ, z being a specific direction in their simulations, at z = 0. Thus, there had been different conceptions in dealing with Neumann boundary condition. However, now, by the method that References [15,45] present, there is a well-motivated mathematical definition of imposing Neumann B.C. as determining the value of Π ab on the boundary. Therefore, there is at last a clear definition of "imposing Neumann B.C." which could guide us through the complexity of GR formalism.
As it is discussed before [28,29,56], assuming Neumann B.C. could lead to three different possibilities. First, by comparing our results with observations, we might realize that the derived value of c 1 is very tiny. Of course, in this case, we will conclude that Neumann B.C. is not reliable, i.e., Dirichlet B.C. triumphs. Second, one might derive different values for c 1 from various observations, i.e., observations based on lensing, rotation curves of galaxies, CMB, etc. report contradictory results for c 1 . Then, we would conclude that our fundamental assumption, i.e., Neumann B.C., is not self-consistent and thus excluded. The last possibility is that the value of the Neumann constants c 1 which is found from the results of different observations is nonzero and compatible at various scales. In this case, the significance of the Neumann B.C. should not be underestimated.

Discussion on Mach's Principle as Boundary Condition
Using the metric in Equation (3), one could derive the field and conservation equations. These are reported in Section 2. Then, one could immediately see from Equation (11) that the scalar fields are not physically independent. Here, one usually concludes that, if the anisotropic stress is negligible, then the two scalars are equal. In other words, it is presumed that ∂ i ∂ j (Φ − Ψ) = 0 leads automatically to Φ = Ψ. However, simple investigation of this partial differential equation ( PDE ) reveals that a presumption is implicitly made about the boundary condition. Namely, it is assumed that Φ − Ψ vanishes at the boundary, i.e., Dirichlet B.C on the boundary is presumed. On the other hand, by imposing a new boundary condition, e.g., Neumann B.C., the solution to the field equations could substantially change. In conclusion, unless particular boundary conditions are imposed, the gauge conditions do not uniquely fix the potentials. See Reference [82], page 144, for a discussion.
Wheeler [16], too, argues that Einstein field equations, as a system of PDEs, do not suffice to define a solution. These equations must also be supplemented by well-defined boundary conditions which he recognizes as Mach's principle. Other formulations of Mach's principle are also possible which might be considered as mathematically vague. See References [16][17][18]21] for more details. According to Wheeler, the role of Mach's principle, i.e., boundary condition, is to select the solutions of Einstein's equations based on the physics of the system. This method of selecting solutions is familiar in solving Poisson equations (electrostatics or Newtonian gravitation) which generally goes without a specific name. In fact, as the boundary condition, the Poisson equation is typically supplemented by the statement that the potential decreases as 1/r. Then, one could see that the distribution of matter (or electric charge) uniquely determines the general form of the potential. If, on the other hand, another boundary condition is imposed, the behaviour of the potential could be entirely different [83]. In the case of general relativity too, the boundary condition is used to select the physically allowable solutions. For example, when one seeks a stellar interior solution, e.g., Tolman-Oppenheimer-Volkoff solution, one might insist that, as the boundary condition, the pressure goes to zero and the solution matches the exterior solution, i.e., the Schwarzschild metric, at the boundary of the star. In the present work, we implement Neumann B.C. on cosmological perturbation equations to pursue a modified model which is more in accord with the observable universe and, thus, free of dark matter and dark energy. As discussed in Section 5.3, by using Neumann B.C., the homogeneity and isotropy of the background universe is still preserved.
However, we do not consider the details here, we only mention that, by imposing Neumann B.C., the super-horizon solutions would also be modified. Moreover, if the parameter c 1 is time dependent, then the results could be substantially different from that of References [7,84] depending on the rate of change in c 1 . In this case, the Weinberg theorem should be modified. See Reference [85] for other models which violate the Weinberg theorem.

Discussion on the Averaging Procedure
We should note that the averaging procedure defined in Section 3 is carried out on a hypersurface of constant time t. In other words, it is presumed that the averaging integrals are performed instantaneously on the whole of the universe. See Amendola and Tsujikawa [11], page 294, for a discussion on this matter. Of course, one could assume that, in the absence of walls and other huge inhomogeneities, the structures in the cosmos form quite similarly everywhere; thus, the presumption of an instantaneous average is somehow justified. However, as an alternative, the averaging procedure could be done by calculating integrals over the "light-cone" which is a more subtle method. See Reference [86] for a complete description of this procedure. The integration in this method is carried out on a section of spacetime which is "causally" connected with the observer (us). Basically, this type of averaging could lead to quite different consequences (at least for some redshifts and some particular mass distributions) compared to the procedure that we used in Section 3; however, a survey on the difference between the two approaches is beyond the scopes of the present work.
The problem of averaging perturbations at large distances needs to be examined more carefully. In this regard, one might argue that either (i) the perturbed metric, when averaged out on large scales, must result in an isotropic and homogeneous background or (ii) the spacial average of the perturbations needs to vanish on large scales. This latter condition presumes that the perturbed metric provides no contribution to the background metric of the universe at large scale. From an observational point of view, the condition i must always be satisfied. We will call (i) and (ii) the weak and strong smoothness conditions, respectively, and it is the intention of this part to show that a constant c 1 satisfies the strong smoothness condition while a time-dependent Neumann parameter c 1 (t) is compatible with weak smoothness condition.
To see this, one could readily check that, by using the definition of A in Section 3, the averaged metric could be found as ds 2 = −dt 2 + R 2 (t) (1 + 2c 1 (t)) δ ij dx i dx j (44) in which we have recovered the time dependence of c 1 for the sake of completeness of the discussion. The elements of this averaged metric only depend on time while the governing cosmic equations, which are reported in Appendix A, are modified now. Therefore, a time-dependent Neumann parameter changes the form of the Friedmann equations by providing a share into the background equations. However, a suitable rescaling of the scale factor as R (t) = R(t) (1 + 2c 1 (t)) 1/2 shows that the geometry of the metric (44) is equivalent to FLRW geometry at any time. In general, we conclude that a time-dependent Neumann parameter satisfies condition i. Moreover, regarding the pure time dependence of the averaged metric in Equation (44) and the assumption that the anisotropic scalar is negligible, i.e., π S = 0, one concludes that the background universe is essentially a perfect fluid.
On the other hand, if c 1 is a constant, then the form of the Friedmann equations remains intact as reported in Equation (27). In this case, which is the main focus of this work, the strong smoothness condition is satisfied. Of course, the result in this latter case is again a perfect fluid. The appearance of the new term, i.e., c 1 (t), in the metric only changes large-scale systems. The reason is that we have only changed the boundary condition of the gravitational field. In particular, in the same way that the expansion of the universe does not change the laws governing the nongravitational forces, i.e., electromagnetic, etc., these laws remain intact under appearance of c 1 (t) in Equation (44). For example, a typical scattering between subatomic particles happens in a very short time compared to the rate which R(t) or c 1 (t) change significantly. Therefore, in calculating a cross section, it is safe to neglect the time evolution of the universe.

Conclusions
In this study, we built a model based on the idea that local scale physics might be affected by the global expansion of the universe through a term which is related to the de Sitter scale of acceleration cH 0 . To do so, we used Neumann B.C. instead of the usual Dirichlet B.C to solve Einstein perturbation equations. This new boundary condition is mathematically well motivated, however, its mathematical well posedness is yet to be surveyed.
The outcome of the model is reminiscent of Mach's principle. The Mach's principle argues that the condition of the "distant objects of the cosmos" somehow enters into the laws of local mechanics. There are various interpretations and formulations of this principle which are thoroughly presented in Reference [18]. See also References [21] or [87], pages 543-549. However, the present model is based on the idea that "the expansion of the universe" is essential and that the new term enters the equation of motion through presuming Neumann B.C. in this expanding universe (Wheeler's formulation of Mach's principle [16]). By imposing a boundary condition, according to Wheeler, one essentially "selects" a specific solution. The Neumann boundary condition which is used here presumes a nonnegligible effect of the distant objects on local dynamics. Another possible path to the present model might be through evaluating the surface term of the generally covariant integral formulation of GR [22] by imposing a Neumann boundary condition.
In the standard model of cosmology, the accelerating expansion of the universe is provided by assuming a constant density of energy due to Λ. This term provides a repulsive force at large scales in ΛCDM while presuming that a dark halo provides a more attractive force (of course through Newtonian force of gravity) at galactic scales. On the other hand, the new term in the present theory, i.e., the term proportional to c 1 cH(t) in the Poisson equation (Equation (20)), provides more attractive forces at galactic scales while still produces a negative deceleration parameter at later times, i.e., see q < 0 in Figure 4. Also, the total equation of state parameter w t tends to w t < −1 at later times, as was shown in Figure 5. The key to understanding this seemingly odd behavior, i.e., making attraction and repulsion at different scales, is that, unlike the ΛCDM model with its static Λ, the new term in our model changes with time because it is proportional to Hubble parameter. As is shown by References [47,48] and here, the existence of a term proportional to H(t) in Friedmann equations could be considered as a model of dark energy.
Unified models of dark matter and dark energy, which rely on "dark fields" instead of dark matter or Λ term, have been proposed before. See [11], chapter 8 and references therein, for a brief review on these models. Among these models, one could mention generalized Chaplygin gas model [49], k-essence model [50], and models based on Bose-Einstein condensation [51]. However, these models might encounter problems in dealing with observations at early or late times. For example, the sound speed in a generalized Chaplygin gas model is small at early cosmic times while it shows a growth at later epochs. Compared with the behavior of c 2 s in the MOD model and ΛCDM as reported in Figure 6, this evolution of sound speed in generalized Chaplygin gas results in incompatibilities with observations of large-scale structure [11]. However, the observational success of the standard model seems to mostly rely on the fact that ΛCDM needs only a few parameters to fit the cosmological observations [9]. Beside the parameters which are usually derived by fitting cosmological data to ΛCDM predictions; there are some other parameters presumed at galactic scales. For instance, a NFW profile (introduced by Navarro, Frenk, and White [88]) relies on two free parameters, i.e., the scale length of the halo and its central density. Furthermore, one also needs more parameters to build a successful particle theory of dark matter, e.g., supersymmetry ( SUSY ), which we will not discuss here. In this respect, the present model could be considered quite economical because at galactic scales it relies only on c 1 . See Equation (20). At cosmic scales too, the Neumann model is dependent to c 1 , ε and Ω c 1 , which seems quite parsimonious. Also, it attributes the puzzling phenomena of dark matter and dark energy to the boundary conditions of the field equations (the state of matter distribution at the boundary) or to the presence of new substances or other dark fields.
Admittedly, the question of boundary/initial conditions in cosmology has not been thoroughly studied so far. The same is true about a possible link between local and global dynamics. However, some exceptions could be found in the works of the founders of the standard model. For example, Dicke and Peebles [89] criticize the proposal of Pachner [90,91], who suggested the existence of a connection between local and global dynamics based on Mach's principle. In fact, Dicke and Peebles argue that the apparent connection proposed by Pachner is only formal, i.e., such effects would be unobservable. In addition, Peebles [1] (part III, pages 361-363) argues that the notion of initial condition is not necessary for a cosmological model unless probably for the fine-tuning problem of very early universe. However, these works do not include a systematic imposition of a new boundary condition and its local effects and they mostly argue, on general grounds, that such modifications would be unnecessary or unobservable.
The degree of reliability of Neumann model needs to be further studied in galactic, cluster, and cosmic scales. At galactic scales, the local and global stability of the systems have to be surveyed. Specifically, by developing the work of Reference [56], one could derive the local stability of a gas + stellar system and compare the results with star formation rate. Moreover, the scaling rules of the MOD model at galactic scales show interesting properties which will be reported in the future. Another challenge is the problem of dwarf galaxies which is well established within the MOND paradigm but seems quite problematic in ΛCDM [36].
However, the main issue at cosmic scales remains to be the CMB analysis of MOD. Regarding the analysis which we presented in Section 2, one could see that the MOD model works based on a modified Poisson equation which is indeed of fourth order. See Equation (3) of Reference [56] for the derivation. The homogeneous form of this equation is usually discussed in the theory of linear elasticity and has been named biharmonic equation. The solutions to biharmonic equation include the solutions of Laplace equation. By implementing the solutions of the modified Poisson equation into perturbation equations, one can derive CMB spectrum. This is a delicate issue which would be reported in a following work. However, as the authors of Reference [55] have discussed, there is already a powerful technique to systematically search for a deviation between the two scalar potentials at cosmic scales, i.e., the so-called anisotropic stress function η = Ψ/Φ, which hints to possible need for a modification in gravity (at 95% confidence level). Acknowledgments: H.S. is grateful to Pavel Kroupa for his support during a visit at Bonn University and for also suggesting some tests to confront the present model with. H.S. is also thankful to Vasanth B. Subramani, Hossein Afsharnia, and Mehdi Ebrahimi for various discussions. The authors also appreciate help from David Merritt. It is a pleasure to thank Mahmood Roshan and Neda Ghafourian, who kindly reviewed the draft and provided insightful comments. Many thanks to Yoshiaki Sofue for providing his galactic data freely. We also acknowledge the anonymous reviewer whose comments helped us to clarify the presentation of the work. The appendix on trajectory of massive particles is added due to reviewer's suggestion. This research has made use of NASA Astrophysics Data System.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Time-Dependent Neumann Parameter
Similar to the method that was used in Section 2, one can prove that, by assuming a timedependent Neumann parameter c 1 (t), the next two pure time-dependent equations could be derived from Einstein field equations: while, by using energy-momentum conservation equation, it is possible to achievė As one could see immediately, these three equations are independent, i.e., the third one could not be derived from the first two equation unless one assumes thatċ 1 = 0. This is the key reason for assuming a constant c 1 in Section 2.
It is worthy to note that the conservation equation could be derived "order by order" from GR field equations. However, the key point here is that Equation (A2) is independent from Equations (A1) because we have mixed zeroth-and first-order terms in these equations to obtain all pure time-dependent terms.
On the other hand, if one drops the assumption that there is necessarily two independent equations, then Equations (A1) and (A2) are still mathematically self-consistent. The reason is that one could derive the evolution of radiation (with P r = 1/3) and baryonic matter (with P m = 0) from Equation (A2) as follows: One can then put the results in Equation (A1) to obtain a system of two independent equations of two unknown functions R(t) and c 1 (t). This system could be solved numerically by presuming suitable initial conditions, i.e., determining a(0),ȧ(0), c 1 (0), andċ 1 (0). The results are very similar to the present model at small z as one expects. However, the dynamical system of such a model is quite complicated and we will not consider it here.

Appendix B. Discussion on the Trajectory of Massive Particles, GR Measurement Process, and the EPS Theorem
Light signals, alongside free falling particles, form the basis of measurement process in general theory of relativity [46,59,92]. Therefore, if the trajectory of the light is modified by Equation (19), the measured values for physical quantities would change too. It is the aim of this appendix to show how the change in the trajectory of photons would lead to a modification in the trajectory of massive particles. However, before we analyse this issue, we need to review the measurement process of spacetime intervals, which lies in the center of a constructive-axiomatic approach to GR [46,93].
The accepted procedure of spacetime measurement in GR is known as the method of geometrodynamic clocks [59]. This method had been introduced by References [92,94], and later modified slightly by Reference [95] to eradicate the need for atomic standards which tend to change when placed in gravitational fields. See Reference [59], chapter 5, for a thorough review on measuring spacetime intervals in curved spaces by using geometrodynamic clocks. These clocks are made of two mirrors with a constant separation, with light rays bouncing back and forth between them while both mirrors are in free fall. See Figure A1 for a sketch of this clock. In a gravitational field, all particles fall at the same rate; thus, geometrodynamic clocks are in fact matter independent. Although the tidal forces might destroy the parallelism of the mirrors, one might continually re-parallel them by using extra clocks.
The key advantage of the method of geometrodynamic clocks lies in using free falling particles and light rays. This is crucial because as Ehlers, Pirani, and Schild have proved [46], the light rays build a conformal structure while free falling particles determine a projective structure. This is known as the EPS theorem, which is inspired by previous works of Hermann Weyl on the foundations of differential geometry and their relation to physics [96]. The EPS theorem could be used to derive many interesting results. Most notably, by applying geometrodynamic clocks, one could obtain affine and metric structures of spacetime manifold. To prove this, the authors of Reference [46] have assumed some well-motivated postulates from which one could define exact operations to measure length and time intervals. See Reference [46] for the exact definitions of conformal, projective, affine, and metric structures. In addition, Ehlers [93] provides details of a constructive-axiomatic approach to GR which employs EPS theorem as the cornerstone of the theory. Einstein field equations are derived in this axiomatic method from a few well-motivated set of postulates. Thus, the EPS theorem is essential in axiomatic description of GR.
The notion of geometrodynamic clocks is also important from a physical point of view because it can be proved that the particles in free fall move along the geodesics of a metric g µν which is measured by geometrodynamic clocks. See Reference [59], page 203, for a proof of this theorem. Therefore, if due to the change in boundary conditions the governing equation of motion of light rays is modified-as it has happened now according to Equation (19)-then one could expect that the measured values of time, space, acceleration, etc. would change accordingly.
In this appendix, we resume the velocity of light c in our equations to compare the orders of magnitudes more easily. Assume a massive test particle P which moves from point A to point B, with an infinitesimal distance d x p and within infinitesimal time dt p . See Figure A1. The proper time for this motion could be written as The particle is accompanied by a geometrodynamic clock, which is meant to measure dt p and d x p . The unit of time t c is considered to be the interval between two ticks of the clock while the unit of length x c is the distance that light goes back and forth in one unit of time, i.e., x c is twice the distance between the two mirrors in Figure A1. The clock measures the time interval dt p and the distance d x p based on these units as explained by References [59,95]. For the light signals of the clock, we need not to worry about the deflection from the radial direction because the background metric is homogeneous and isotropic; therefore, any deflection would have to be of first order and its effect in the perturbation metric would be of the second order. Thus, any first order change in the geodesic of the light signals would be toward the center of the perturbation. For any light ray within the horizon (including the signal of the clock) which is moving toward the perturbation, one could write in which x r c stands for the infinitesimal displacement in the radial direction of the signal. We will use this equation in the following to derive the effect of imposing Neumann boundary conditions on the clocks. A B P Figure A1. Trajectory of a particle P in the field of a massive object is measured by using geometrodynamic clocks. The light rays of the clock are shown by solid lines. The clock itself is affected by the gravitational field of the central mass and the boundary condition. Accordingly, the trajectory of the massive test particle too will be modified because the physical quantities such as velocity, acceleration, etc., are measured by geometrodynamic clocks.

Appendix B.1. Trajectory of Massive Particles
In this part, we first prove that imposing a Neumann boundary condition on perturbation equations would result in negligible terms on the right-hand side of the geodesic equation, i.e., −Γ i µν q µ q ν , for massive particles. Then, we will show that, because imposing Neumann B. C. would affect the light signals of the geometrodynamic clock, the left-hand side of the geodesic equation of massive particles dq i /dτ will ultimately change to include a term proportional to cH. Also, we will provide a simple observational evidence for the existence of a term proportional to cH in the equation of motion of massive particle. In addition, a proof of Equation (20) is provided at the end of the results. During our calculations, we will use Table A1, which provides the order of magnitude of the perturbation terms, to systematically neglect the terms with insignificant share in the equations.
The new term, i.e., 2c 1 cH, has been introduced into the equation of motion of massive particles by Reference [28] on the general ground that, because the light rays are disturbed by Neumann B. C., its effect must show itself into the equation of motion of massive objects. However, Reference [28] provides no conclusive method of how one could include the change in geometrodynamic clocks (due to Neumann boundary condition) into the equation of motion. However, in that work, a sample of 101 galactic systems has been surveyed to show that, at the last data points, in which the Newtonian attraction is small, the ratio of the observed acceleration to cH 0 is almost constant. In addition, MOND literatures have already revealed ample signs of the reliability of cH 0 to the galactic data [32][33][34][35]. In this sense, the inclusion of 2c 1 cH 0 into the equation of motion had been justified through a phenomenological approach by Reference [28]. Table A1. Typical magnitudes of perturbations in solar, galactic, and cluster scales. We have assumed typical rotation velocity of the objects in the solar system (galaxy and cluster of galaxies, respectively) as 30 km/s (100 km/s and 1000 km/s, respectively), mass of the Sun M = 1.988 × 10 30 kg (10 12 M , 10 14 M ), and typical distance to the Sun AU = 1.496 × 10 11 m (10 kpc, 10 Mpc). In general, this table shows that the terms proportional to v/c,Φ, and HΦ are negligible in Equation (A9). Also, solar system is governed by the term ∇Φ (the effect of 2c 1 cH due to boundary condition is very small in this scale [29]) while the effect of 2c 1 cH is considerable at galactic and cluster scales. In the last row, we have reported typical observed accelerations a obs in these scales.

Solar System Galaxies Cluster of Galaxies
a obs a 0 ≈ 10 +8 0.01 to 10 0.001 to 0.1 To find the equation of motion of a massive particle with mass m, we write the physical momentum as while the co-moving momentum q, to first order of approximation, could be found as (Dodelson [97], page 102) Please note that Φ and Ψ in Dodelson's book corresponds to −Ψ and Φ in our notation. Then, the geodesic equations for a massive particle could be found as (in first order of approximation) by applying the Christoffel symbols of Equation (4). Now, by imposing Neumann B.C. Ψ = Φ − c 1 into the last two equations, one could arrive at Here, we need to evaluate the terms on the right-hand side of Equations (A8) and (A9). First of all, we assume particles with non-relativistic velocities | p | /mc 1. Thus, we would have E mc 2 . Now, by considering typical values of perturbations from Table A1, one could readily see that the terms on the rhs of Equation (A8) are negligible. Thus, we could find dq 0 /dτ p 0, which means that E is conserved. The right-hand side of Equation (A9) too could be estimated in the same way. Comparing these terms with the values of Table A1, one could see that the only term which survives on the rhs is the first term, i.e., ∝ ∇Φ. Thus, Equation (A9) could be approximated as Therefore, we find out that the imposition of Neumann boundary condition does not provide a significant share on the rhs of Equation (A9). The same result had been reported by Reference [28]. Now, we evaluate the left-hand side of Equation (A9) to derive the equation of motion of massive particles. Knowing that the proper time of the particle could be written as follows , the left-hand side of the geodesic equation could be found as Since v/c and Φ are considered to be small according to Table A1, the last equation could be simplified as d dt p Checking the values of Table A1 again, one could immediately see that the first term on the left-hand side is negligible; thus, we have in which x p is the physical coordinate of the particle. If we presume Dirichlet B. C., i.e., Φ = Ψ, the above equation would readily provide the equation of motion in Newtonian limit. Now, we want to prove that, because the trajectory of light signals is altered, the effect would adjust acceleration measurements in geometrodynamic clocks. Suppose that a geometrodynamic clock accompanies a particle, as in Figure A1, on a trajectory which is parametrized by λ. The time t p and displacement x p of the particle are measured based on the units of the geometrodynamic clock. The clock measures time intervals directly t P = Nt c , i.e., the motion has occurred in N units of time t c . However, the displacement will be a function of both t c and x c , i.e., x p = x p (t c , x c ) because the units of geometrodynamic clocks are always linked by being null, i.e., Equation (A5). We will write all derivatives with respect to λ. Also, as mentioned above, since the deflection from the radial direction is negligible for light signals, we will only consider the displacement of light rays toward center x r c . and while the second term on the second row could be found as Putting Equations (A20)-(A22) together, one could find in which ,in the fifth row, we have imposed Neumann boundary condition Ψ = Φ − c 1 while, in the last row, we have used the fact that cHΦ is negligible in local scales according to Table A1. Finally, by employing this last equation, the equation of motion could be written as or in physical coordinates x p as For a particle around a galaxy, one expects that the centripetal acceleration due to Newtonian potential decreases while there is still a constant acceleration term due to the boundary condition. Thus, if we measure the acceleration at the outermost points of galactic rotation curves, the Newtonian share would be minimized and we could have an estimation of the parameter c 1 . To do so, we have use data from Reference [98] that contain 551 galaxies (including three subsamples). The result is provided in Figure A2. In this plot, R last shows the radius of the last data point (in kpc) while V last is its rotational velocity (in km/s). We have assumed H 0 = 70 km/s/Mpc. For galaxies with small R last (these are the systems for which the observations of rotational velocity are not extended enough), the acceleration due to Newtonian potential is still dominant. Thus, the ratio of (V 2 last /R last )/cH 0 is high for these systems. The maximum value of 2c 1 is found to be about 1. As we move toward right of the plot-systems with extended rotation velocities-the value of 2c 1 becomes more concentrated in the interval 0.01 < 2c 1 < 0.13 (88 % of the data points). Eight systems (about 1.5 % of the data) show a value less than 0.01 for 2c 1 . Three galaxies out of eight posses very large R last . These eight systems are probably galaxies with strong interaction with their neighbouring environments (at least at their outer parts). See also Tables 3-5 of Reference [28] for a similar analysis of 101 galaxies. The value of the Neumann parameter is therefore estimated to be with the bound of 0.01 < 2c 1 < 0.13. This coincides with some previous estimations [28,29,56] Figure A2. Estimating 2c 1 based on the last data points of the rotation velocity. Systems with small R last are more probably governed by the Newtonian term; systems with very large R last might be affected by the long-range interaction of their environment, which is proportional to ∝ 2c 1 cH 0 . On the other hand, galaxies with extended rotation velocities 20 kpc < R last < 60 kpc show a concentration around the value 2c 1 = 0.06. One could check that about 67% of the data points lie within the bound of 0.03 < 2c 1 < 0.12. Data of this plot is extracted from Reference [98] and includes three subsamples (containing 551 objects): S-sample provided by Reference [99], P-sample which contains galaxies with W1-band photometry [100], and C-sample (galaxies with r-band photometry and rotation velocities reported by Courteau [101,102]).
Since, the long-range interaction is always present in the MOD model, it is essential to go beyond the simple approximation of point-like perturbation. To survey systems with many sources of gravitation, one should sum (integrate) over different patches of matter to derive [29] in which x p still represents the physical units. The factor c 2 on the left-hand side makes Φ dimensionless as it is assumed in Section 2. This potential could lead to the next modified Poisson equation: In addition, to write the integrals in the co-moving coordinates, which is used in Section 2, we should note that dx p = R(t)dx and ∂ x p = ∂ x /R(t); thus, the modified Poisson's equation in the co-moving coordinates could be written as which is equal to Equation (20). One should note that, even when the particle is moving away from the perturbation, for which the radial direction of the light rays should be assigned with a + sign on the right-hand side of the Equation (A5), the result for the equation of motion would still be the same. The reason is that, according to Equation (19), the imposition of Neumann boundary condition always leads to more gravitation so far as c 1 > 0. The direction of the light signal does not change this point. In fact, Reference [28] presumes c 1 > 0 to solve the mass discrepancy problem in lensing systems.
Also, we should point out that the effect that a boundary condition introduces to the geometrodynamic clocks could not be removed by readjusting the clocks. The reason is that, unlike the influence of the tidal forces which could be eradicated by re-paralleling of the mirrors, the effect of the boundary condition happens even in a single tick of the device. Thus, it must be included within the equations of motion.

Appendix B.2. The EPS Theorem and the Structure of Spacetime
As we proved in this appendix, because the velocity of a massive particle at galactic scales is typically much smaller than the speed of light, the modification in the geodesy of the massive object (due to Neumann B.C.) is far smaller than cH. However, as discussed above, since the trajectory of the light signal includes a term proportional to de Sitter scale of acceleration, the measured values of the massive particle's acceleration would include the term cH. This procedure is motivated by the observational fact that both massless and massive particles trajectories show dependency to cH. A good example for the former is the observations of strong lenses [28], while for the latter, we could mention the rotation curve data of galaxies [29], local stability of galactic disks [56], and Figure A1 here. In this way, the missing mass discrepancy is united by introducing a new term into the equation of motion.
On the other hand, if one presumes the other procedure of spacetime measurement, e.g., atomic standards [59], then the two theorems of EPS [46] and Ohanian [59] could not be used. Thus, the unification between the massive and massless particle's trajectories would not be possible and the missing mass discrepancy would remain unsolved for massless particles. To summarize, the inclusion of the notion of geometrodynamic clock is the second critical assumption of Neumann cosmology (beside Neumann B.C.).
One could even further argue that the viability of Neumann cosmology based on observations at different scales could be considered as a direct "empirical evidence" for the advantage of the Marzke-Wheeler [92] measurement method (geometrodynamic clocks) over other methods of spacetime measurement. Especially, this could demonstrate the correctness of EPS theorem (within the present model), which is a key tool in constructing the geometrical structure of spacetime.
one might simply argue that these are meaningless numerical coincidences. However, in light of the present work, these dimensional relations might also be interpreted as evidence that the microphysics is related to the cosmological dynamics through Mach's principle.
Knowing that the Hubble parameter in Equation (A31) varies with time, if one assumes that the relation in Equation (A31) is indeed fundamental, then at least one of the other constants in this equation must change with time. The factor G is most commonly speculated as the time dependent parameter in the literature; otherwise, one must be open to reformulating all atomic physics by introducing c, m π , or h p as time-dependent variables. Many interesting results have been derived in the literature based on G(t) strategy. See, for example, Weinberg [105], chapter 16, for an interesting cosmic model based on a time-dependent gravitational "constant". The assumption that G is time dependent is the core idea in Mach1 formulation of Bondi and Samuel [21].
However, the Neumann model presents another possibility since, in all our formulation, the new acceleration term is introduced as c 1 cH 0 in which the parameter c 1 could be, in general, time dependent. The reason is that c 1 represents the effect of the distant objects, which are themselves evolving. Thus, it might be insightful to introduce c 1 (t) in Equation (A31) and to investigate its time variability (instead of G). However, as mentioned before, the Friedmann equations would be modified in this case. See Appendix A for the details. Assuming the Compton wavelength of pion h p /m π c, one might rewrite Equation (A31) as If the left-hand side of the above equation is assumed to be independent of time, then we could see that there should be a relation as c 1 (t) ∝ 1/H(t) to maintain the rhs too as time independent. Therefore, the dimensionless parameter c 1 (t) must be very small in early epochs, i.e., when H(t) is very large, and then show its effects in later times. This provides a new possibility for cosmic evolution; however, for the sake of brevity, we do not follow this case here.

Appendix D. A Semi-Newtonian Model of Gravity: Some Properties and Issues
Having built a modified dynamical model based on imposing Neumann boundary condition on cosmological perturbation equations, a question naturally arises about the possibility of such a model in Newtonian mechanics. At first glance, it seems that there is no such possibility in Newtonian dynamics because this theory exploits only one potential, i.e., there is Φ but no Ψ. Therefore, even in a Newtonian model of cosmic expansion, it seems that one could not derive a term proportional to a 0 = cH 0 due to the limited degrees of freedom. It could be also proved that such an exotic term does not appear in the equation of motion even by imposing Neumann B.C. on Poisson equation. Assuming the classical Poisson's equation with Neumann B.C., in which ∂Φ( r) ∂ n = g( r) is the normal derivative at the boundary of the volume ∂V, one could see that the normal derivative of the potential is related through the divergence theorem to its Laplacian, i.e., where S is the surface element on the boundary. See Reference [106], pages 619-620, for more details. In this section, the potential has the dimension of Length 2 /Time 2 as it is usually presumed in Newtonian mechanics. This last equation is Gauss's law in gravitation, and it imposes a strict solvability conditions on the Poisson equation with Neumann boundary condition. Thus, one cannot choose quite arbitrarily the value of the force on the boundary because it is related to the total mass of the bulk, i.e., |∂Φ( r)/∂n| = GM in /r 2 , in which M in is the total mass enclosed in a sphere of radius r. In other words, no fifth force can appear from imposing Neumann B.C. in Newtonian gravity. For the sake of completeness, we mention that it is always possible to rewrite the Poisson equation as whereρ is the average density of the bulk. This equation is always solvable with Neumann B.C. [106]. However, it should be mentioned that this modified Poisson equation cannot solve the dark matter problem since it decreases the amount of gravitating source byρ. However, it is still possible to construct a semi-Newtonian model of local dynamics similar to the model introduced in Section 2 by considering Mach's principle. To do so, we will use a well-defined property of classical mechanics and three well-supported observational facts regarding our universe. In classical mechanics, (a) Newton's second law is valid only in inertial frames. The observational facts used here are as follows: (b) All astrophysical systems are accelerating; thus, they could only be considered as an approximation of inertial frames. (c) The typical acceleration in the scale of galaxies and cluster of galaxies is of the order of de Sitter scale of acceleration cH 0 . (d) The universe is isotropic and homogeneous on very large scales. Assumptions a, b, and d are well accepted within the community while assumption c is supported by an impressive wealth of observations modeled based on the MOND paradigm. Now, imagine a particle moving around a galaxy. See Figure 1. This galaxy is not, in principle, an inertial frame; however, we could assume an inertial frame i far from the galaxy g and could write the equations of motion of a particle p, with mass m, in that inertial frame. Then, the position of the particle with respect to the inertial frame could be written as r pi = r pg + r gi ; thus, the Newton's second law in that frame is as follows: m( a pg + a gi ) = −m ∇Φ N (A37) in which Φ N must be derived from Poisson equation, | a pg | = v 2 /r provides the centripetal acceleration of the particle p around the galaxy g and a gi is the acceleration of the galaxy with respect to the inertial frame. In writing Equation (A37), we have used the assumptions a and b above. Also, based on assumption c and to solve the mass discrepancy problem at galactic scale, we presume that the value of a gi is of the same order of de Sitter acceleration. This assumption also helps us to understand how far the inertial frame ought to be; however, the argument comes from the causal structure of special relativity and as a result of a time-dependent boosts. In fact, for an accelerating observer with acceleration g, the local frame attached to the observer could be safely used to set coordinates within the Rindler horizon defined as the distance L R ≡ c 2 /g. For example, for typical accelerations on earth, i.e., g ≈ 10 m/s 2 , one has L R ≈ light year which is quite satisfying. However, at galactic scale with a typical de Sitter acceleration cH 0 , the Rindler horizon is of the order of the Hubble distance L R ≡ c/H 0 . Thus, one infers that the origin of the term a gi is related to the general evolution of the cosmos. Since the cosmic fluid at Hubble distance is homogeneous and isotropic, i.e., assumption d, the only form of a gi compatible with this condition is when a gi is radial. Thus, the magnitude of a gi is of the order of cH 0 and it is toward the center. Therefore, the governing equations in this semi-Newtonian model could be written as a pg + 2c 1 cH 0ˆ r = − ∇Φ N (A38) in which the term 2c 1 is introduced to make this model compatible with the Neumann model of Section 2. It should be noted that no horizon, including Rindler horizon, exists in Newtonian mechanics and that the inertial frames of Newtonian physics are presumed to be satisfied quite globally. The notion of Rindler horizon is borrowed from special relativity, hence, the name semi-Newtonian. It is also worthy to note that the result derived in this part fulfils Mach3 formulation of Bondi and Samuel [21], who imply that the "local inertial frames are affected by the cosmic motion and distribution of matter"; however, Mach3 is usually discussed in connection with rotating frames. In conclusion, the inertial frame i is not attached to a "single object" very far away; i is linked to the homogeneous and isotropic matter which is distributed at cosmic scale. In this way, the symmetry of the model is preserved. The appearance of the new term on the left-hand side of Equation (A38) stresses the fact that this model is indeed a modified dynamical model. However, to use the full capacity of potential theory, it is more appropriate to transfer this term to the rhs of the equation of motion and to rewrite Equation (A38) for a system of particles as a pg = − ∇Φ (A39) in which Φ is the modified potential which must be derived from the second equation, i.e., the modified Poisson equation. The method to derive this modified Poisson equation is exactly the same as the method used in the case of Neumann model in Reference [29].
It is important to note that this semi-Newtonian model shares all transformation properties of the Newtonian theory. To see this point, assume the next extended Galilean group of transformations in which R ij and d i are time-dependent rotation matrix and vector, respectively. See Reference [107] and [93] for more details. Then, for the Galiliean group G, we haveṘ ij = 0 andd i = 0 while, for the Newtonian group N , we haveṘ ij = 0, thus, G ⊂ N . As References [93,107] have explained, under a Newtonian transformation, the gravitational potential transforms as Φ( r, t) → Φ( r , t ) = Φ( r, t) − r . A (A41) in which A is some constant acceleration. This transformation explains in an elegant mathematical way that the field strength is frame dependent, i.e., it changes by an amount of − A under a Newtonian transformation (the basis of Einstein's elevator argument). On the other hand, when the matter distribution is given, the tidal forces could be solely determined. The semi-Newtonian potential of Equation (A39) preserves the transformation in Equation (A41) as one may check easily. As a result, the homogeneity of the universe will be preserved under this model. This is a crucial requirement for models of modified dynamics. As Benisty and Guendelman [108] have argued, when Newton's second law is modified by a nonlinear function of the acceleration, as the MOND model proposes, the concept of relative acceleration would be lost. As a result, it is probably not possible to maintain the homogeneity principle anymore.
Although the semi-Newtonian model seems to be a mathematically consistent model, we have clear reasons to favor the general relativistic view of Section 2. First of all, the appearance of the new acceleration term in Equation (A37) had to be "introduced" based on four assumptions while such term is "derived" in Section 2 by imposing Neumann boundary conditions and as a result of background-perturbed metrics coupling. (Of course, to find a modified equation of motion at a large scale which goes beyond Newtonian mechanics, one must be open to introduce some aspects of the phenomenology of large-scale systems into the model [90].) Moreover, even if we attribute 2c 1 cH 0 in Equation (A39) to distant objects (Mach's principle), the notion remains conceptually vague in contrast to Wheeler's idea, which is mathematically well defined. Furthermore, we had to use the the concept of Rindler horizon to determine the domain of reliability of inertial frames in Newtonian mechanics. Regarding this last issue, the basis of the semi-Newtonian model could be strengthened by considering accelerated frames in special relativity. In this way, the Rindler horizon would naturally emerge. See Reference [30] for a relativistic method of introducing de Sitter acceleration into the equation of motion and its consequences at galactic and extra-galactic scales [31].
The second reason to favor the GR version of modified dynamics is the clear triumph of the Einstein equation in explaining many gravitational phenomena specially in the case of strong gravitational fields. This could not be incorporated in a semi-Newtonian model of gravity. Third, precise experiments of high-energy collisions, time dilation effect, spin-statistics theorem, etc. ensure us the superiority of the Lorentz group of special relativity over the Galilean group of Newtonian mechanics. Thus, as a natural extension of special relativity, GR provides a better chance toward unification of physics.