1. Introduction

The universe has been extensively studied by cosmologists, astrophysicists and many other scientists. It is an immensely huge body with observable volume estimated to be in the order of 4 × 10⁸⁰ m³. Although there are some 2 billion galaxies with as many as an estimated 1 × 10²⁴ stars, the universe is only sparsely populated by matters with an estimated critical density of 10⁻²⁶ kg/m³ (or 10 hydrogen atoms per cubic meter) [1]. With such a low average density it is easy to see the difficulties of developing an effective model for the universe. In additional, known material objects vary from sub-atomic particles to stars with size of thousands of solar radius, presenting a mathematical challenge especially in numerical analysis. Many believe that the space between those known heavenly objects are filled with particles that have zero mass but with momenta governed by various not so well defined fields [1]. In more recent time, it has been suggested that the vast empty space could be filled with yet unknown dark matter and dark energy [2]. If this is the case, it is even more difficult to develop a workable model as most of the properties of dark matter and dark energy are unknown.

In this paper, we intend to bypass all the complexities of how to set up a set of field equations for the universe. Instead, we use the fact that we could use information given to us through electromagnetic waves we are receiving from the universe all the time. Various properties of electromagnetic waves have already been used to measure or estimate distance, movement, mass or chemical compositions of distant bodies in the universe. Through measurement of cosmic microwave background radiation, the background temperature of the universe has been determined to four-digit accuracy [3]. We intend to use this information to set up a field equation so that the temperature histories of the universe could be predicted.

With the assumption of the cosmological principle that the universe is homogeneous and isotropic everywhere, the temperature could be found from a single energy balance equation in that the energy loss through electromagnetic waves is assumed to be the product of an energy transfer coefficient and the temperature itself. Other forms of energy production or consumption are altogether assumed to be represented by a single energy source. For examples, nuclear fusion could produce energy while gravity involves negative energy. Applying the cosmological principle, all the processes involved anywhere in the universe could be represented by a single average term. We shall show how parameters could be fine-tuned to give various temperature histories.

Our research is motivated by our belief that the cosmic microwave background (CMB) needs not come from a very high temperature at the distant past. The predominant current thinking is that this radiation began when the universe had cooled down substantially from its initial infinite temperature but still dense enough (about the density of air) to have sufficient number of hydrogen atoms (or other form of matter) to scatter photons. These photons are what we are mapping now. The models we used in this paper do not follow the same line of postulation. We consider that CMB is that of a black body radiation (also known as thermal blackbody emission). This process has been the same throughout the whole history of the universe. We also measure the amount of energy emitted as per unit mass. It is, therefore, not directly dependent on density. Despite the fact that our models are so simple in terms of mathematics, we show that our models have the flexibility to cater for a wide range of demands including the initial temperature and net energy source. We also show how the nonlinear nature of CMB could be managed numerically without involving the uncertain process of backward marching in time.

Although it is not related directly to how the universe could be modelled, we have chosen a solution of the nonlinear Schrödinger equation to demonstrate how optical vortex solitons could be generated if conditions are suitably chosen [4]. Vortex solitons have been observed experimentally [5] [6]. A short explanation of the Schrödinger equation is that it represents equilibrium between energy loss through the linear dispersion terms and energy gain through the nonlinear terms. It is unlikely that electromagnetic waves will travel through space in the same form because in an optical fibre the waves are confined within the fibre itself. But we are dealing with the same electromagnetic waves. The difference is that, in the space, we may not find the set of conditions that is needed for this particular mode of transmission. Nevertheless, we believe that this example will provide us with useful insights on how mathematical modelling could be used.

2. An Energy Balance Equation

Consider that the time rate of generation of all forms of energies per unit mass can be represented by a single term, Q in J/kg/s. Disrespecting to what form or forms by that energy are being removed, the rate of total energy loss per unit mass for the system is assumed to be proportional to the system temperature, T. Then, the rate of temperature change is given by the following equation:

$c_m\frac{\partial T\left(t\right)}{\partial t}=-hT\left(t\right)+Q$ (1)

where c_m is the specific heat capacity in J/K/kg and h is the energy transfer coefficient in J/kg/K/s. It should be noted that written in this form, density is not involved and h is a constant characteristic of the system as a whole.

Using scaling factors, S_T in Kevin for T and S_t in second for t, both temperature and time could be reduced to dimensionless quantities, T^* and t^*:

$T^{*}=\frac{T}{S_T};t^{*}=\frac{t}{S_t}$ (2)

Then Equation (1) could be rewritten into a dimensionless form:

$\frac{\partial T\left(t\right)}{\partial t}=-c_{2}T\left(t\right)+c_1$ (3)

where the superscript * for both T^* and t^* has been omitted for simplicity reason, and

$c_1=\frac{Q{S}_t}{c_m{S}_T};c_2=\frac{h{S}_t}{c_m}$ (4)

3. Models Based on Linear Assumptions

For linear models, c₁ and c₂ in Equation (3) are assumed to be constants. Then the solution of Equation (3) is

$T\left(t\right)=\frac{c_1}{c_2}+\left[T_0-\frac{c_1}{c_2}\right]\exp \left(-c_{2}t\right)$ (5)

where T₀ is the initial value of T. As it is not possible to take measurement of the universe as a whole and the various physical variables used in Equation (1) are unavailable, it is assumed that c₁ and c₂ in Equation (4) could be determined by calibration or fine tuning as illustrated in the following models:

1) A stationary model

In this model it is assumed that the universe has reached its stationary temperature so that,

$T\left(t\right)=\frac{c_1}{c_2}$ (6)

That is, at long enough time, the second term, the contribution of the initial value in the solution, Equation (5), is no longer significant. From observing the cosmic microwave background, this background temperature of the universe has been measured at 2.73 K [2]. This information could be used to fine-tune the system. Choosing S_T = 1 K and converting to dimensional quantities, then from Equation (4),

$\therefore \frac{c_1}{c_2}=\frac{Q}{h}=2.73\text{K}$ (7)

As c₁ and c₂ are found only as a ratio, we still cannot use Equation (5) to work out the histories as well as the initial temperature. Other assumptions are needed as described in the next model.

2) A zero net energy model

This popular model is based on the argument that although the universe contains billions of energy generation bodies such as our sun, the gravitational forces needed to keep the system together with other energy consuming processes use up equal amount of energy resulting in zero nett energy production. The reason that the universe still has a small temperature of 2.73 K is due to the remnant of its extremely high temperature past.

For a zero net energy universe, the constant, c₁ = 0, in the field Equation (3) and also in its solution Equation (5). But, since T₀ the initial value could be chosen arbitrarily, c₂ can only be found by knowing a particular temperature at a specific time. The widely accepted assumption, based on distances of faraway stars or on the rate of expansion of the universe, is that the universe has an age of 13.8 billion years. There are widely different suggestions about the initial temperature; as an illustration, it is assumed that T₀ = 10¹⁴ K. Based on S_T = 1 K and S_t = 1 year and using Equation (5), c₂ could be found by solving:

${10}^{14}\times \exp \left(-c_2\times 13.8\times {10}^9\right)=2.73$ (8)

to give c₂ = 2.273 × 10⁻⁹. With c₁ and c₂ known, Equation (5) could be used to give the whole temperature history.

3) A composite model

The present observed temperature could be assumed to have come from both the heat source and the remnant of past temperature. In all cases, the system parameters, c₁ and c₂ could be fine-tuned so that the universe could start from any given initial value and cool down or heat up to the present temperature after a specific period of time. If T₀ = 10¹⁴ K and the age is 13.8 billion years and the presently observed temperature of 2.73 K is coming from a combination of the remnant of T₀ and heat source Q at different proportions, the system needs to be fine-tuned according to the values given in Table 1, where Percentage refers to the percentage contribution from the remnant of T₀ and the last two columns are what the system parameters should be fine-tuned to give the present observed background temperature of 2.73 K.

As the time span involved is so large, it is not a surprise to see that c₂ is varying over such a very narrow range. In fact, it is only from analytical solutions that c₂ could be so determined.

4) Non-linear Models

Although the linear models could be fine-tuned to fit what have been hypothesized as the temperatures in the universe, namely the background temperature and the age, there is no justification that the system parameters, c₁/c₂ and c₂ could both be constant over such a huge temperature range. Non-linear models using temperature dependent parameters would be more likely to give better predictions. But it is unlikely that analytical solutions could be found for such nonlinear models. We shall use instead a linearized approximation by subdividing the time span into segments so that over each time interval the system parameters are replaced by constants that represent the average over that particular time interval. To fine-tune such a system, it is necessary that the temperatures at the beginning and at the end of each segment are known or could be estimated.

An example with c₁ = 0 (that is assuming zero energy source.) and T₀ = 10¹⁴ at t = 0 and temperatures at each other period are as given, the values for c₂ found are shown in Table 2. It should be noted that the chosen time and temperature values in Table 2 are compatible with those often used in Big Bang theory [7].

These results show that the cooling rate would be very high at the beginning and slowed down to a small value only toward the end of the current time history of the universe. To lower c₂ at the beginning, a plausible model may assume that there is a rapid expansion that consumes so much energy leading to a negative c₁. For example, if c₁/c₂ = −10¹¹, then for the first period between 0 and 180 seconds, the value for c₂ = 1.21 × 10⁶. However, the above is only one of unlimited number of possibilities. In fact, the system could start with a zero temperature providing that the heat source is positive as in the cases shown in Table 3, for c₂ = 1.

5) Low temperature models

The problem with high temperature models is that it is difficult to explain how the temperature could be cooled down so rapidly from infinite to say 10¹⁴ K in just 10⁻⁴² second. In fact, if the initial temperature is infinite as postulated, there is no way that this temperature could be lowered, no matter how fast and how much energy, that is any amount less than infinite, is removed, by expansion or any other means. Because it is difficult to justify both physically and mathematically, dark matter and dark energy, both unknown to us, have been suggested as the best plausible explanations.

It is much easier to assume a low temperature model. For example, taking the energy source to be the same as in Case (a), that is c₁/c₂ = 2.73, and cooling coefficient in Case (b), that is c₂ = 2.273 × 10⁻⁹ the time requires for temperature to

cool down from 1000 K to 3 K is only 114.6 years. If starting at zero temperature, the time taken for the temperature to heat up to 2.70 K is 62.9 years.

It is better still to assume that the starting temperature is 2.73 K, the same as the observed temperature now and Q = 2.73 h, then the universe is always maintained at a constant temperature of 2.73 K all the time.

An often used justification for a high temperature start is based on the second law of thermodynamics. As the temperature is so high, the initial entropy is very low. Then, all spontaneous processes supposed to take place subsequently should bring down the temperature and hence an increase of entropy. But it should be noted that out of all the cases studied above from (a) to (e), with the case of a zero starting temperature as an exception, energy is coming from a higher temperature than the loss term associated with h. The implication has been that entropy is always increasing as requires by the second law.

4. Electromagnetic Waves Transmission

Consider special electromagnetic waves called vortex solitons. Using cylindrical coordinates that is natural for this type of wave transmission, the nonlinear Schrödinger equations [4] in the spatiotemporal domain are:

$\begin{array}{l}i\frac{\partial u}{\partial z}+\frac{1}{2}\left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}\right)+u\ast v-\sigma \left({\left|u\right|}^2/4+2{\left|v\right|}^2\right)u=0\\ 2i\frac{\partial v}{\partial z}+\frac{1}{2}\left(\frac{{\partial }^{2}v}{\partial r^2}+\frac{1}{r}\frac{\partial v}{\partial \theta }+\frac{1}{r^2}\frac{{\partial }^{2}v}{\partial {\theta }^2}\right)-qv+u^2/2-\sigma \left(4{\left|v\right|}^2+2{\left|u\right|}^2\right)v=0\end{array}$ (11)

where u and v are the amplitudes of the coupled electromagnetic waves, q is a coupling coefficient, σ is the nonlinear coefficient and z the distance travelled.

The above equation could be solved by the Lanczos-Chebyshev pseudospectral method [8], seeking solutions in the form:

$\begin{array}{l}u\left(r,\theta ,z\right)={\displaystyle {\sum }_{i=0,M}{\displaystyle {\sum }_{j=0,N}{u}_{ij}\left(z\right)r^i{\theta }^j}}\\ v\left(r,\theta ,z\right)={\displaystyle {\sum }_{i=0,M}{\displaystyle {\sum }_{j=1,N}{v}_{ij}\left(z\right)r^i{\theta }^j}}\end{array}$ (12)

For the given computational 2D domain, an affine transformation is carried out to rescale it into new coordinates $\left\{r,\theta \right\}$ to $\left[-1,+1\right]\times \left[-1,+1\right]$. Spatial-discretization of Equation (11) is carried out by using collocations at points corresponding to the roots of Chebyshev polynomials:

$\begin{array}{l}{r}_i=-\cos \left[\frac{\left(2i+1\right)\pi }{2\left(M-2\right)}\right],i=0,\cdots ,M-2\\ {\theta }_j=-\cos \left[\frac{\left(2j+1\right)\pi }{2\left(N-2\right)}\right],j=0,\cdots ,N-2\end{array}$ (13)

This approach would be the most economical choice numerical-wise; more details about the Lanczos-Chebyshev psudospectral method may be found in Reference [9].

Equation (11), together with the boundary conditions, is then transformed to a set of first order ordinary differential equations in matrix form of

$iA{U}_z-LU-H\left(U,z\right)=0$ (14)

where $U$ is a vector ${\left(u_{11},u_{12},\cdots ,u_{1M},\cdots ,u_{N1},u_{N2},\cdots ,u_{NM},v_{11},v_{12},\cdots ,v_{1M},\cdots ,v_{N1},v_{N2},\cdots ,u_{NM}\right)}^{\prime }$, $A$ and $L$ are linear matrices, $H$ is a nonlinear vector.

Using the unconditional, stable and implicit Crank-Nicholson algorithm, we set

$iA\left(U^{m+1}-U^m\right)-\frac{\Delta z}{2}\left[L\left(U^{m+1}+U^m\right)+H\left(U^{m+1},z^{m+1}\right)+H\left(U^m,z^m\right)\right]=0$ (15)

where the superscript m refers to the time-step number. Since Equation (14) is nonlinear, we let

$U^{m+1,0}=U^m$, then $U^{m+1,r-1}\to U^{m+1,r}$ (16)

where the superscript r refers to the iteration number and the arrow implies that integration is used. It should be noted that the iteration involves only one nonlinear $H$ term and since $A$ and $L$ are linear, only one inversion is required for all the time steps:

$U^{m+1,r}={\left[iA-\frac{\Delta z}{2}L\right]}^{-1}\left\{\left[iA+\frac{\Delta z}{2}L\right]U^m+\frac{\Delta z}{2}\left[H\left(U^{m+1},z^{m+1}\right)+H\left(U^m,z^m\right)\right]\right\}$ (17)

The system is a nonlinear initial value problem. In theory any set of arbitrarily chosen initial values could be used as non-equilibrium components present initially would disappear through the marching process. But some unwanted components may grow instead. Therefore, numerically, to ensure success, it is necessary to select as good as possible a set of initial values and also to use a small enough step size. In our case, we took advantage that an approximate solution could be obtained more easily. If we consider that

$u\left(r,\theta ,z\right)=\exp \left(iS\theta \right)f\left(r,z\right),v=\exp \left(2iS\theta \right)g\left(r,z\right)$ (18)

where $S=0,1,2,\cdots$ is the integer vorticity. Substituting Equation (18) into Equation (11) gives,

$\begin{array}{l}i\frac{\partial f}{\partial z}+\frac{1}{2}\left(\frac{{\partial }^{2}f}{\partial r^2}+\frac{1}{r}\frac{\partial f}{\partial r}-\frac{S^{2}f}{r^2}\right)+f\ast g-\sigma \left(\left|f\right|/4+2{\left|g\right|}^2\right)f=0\\ 2i\frac{\partial g}{\partial z}+\frac{1}{2}\left(\frac{{\partial }^{2}g}{\partial r^2}-\frac{2iSg}{r}+\frac{4S^{2}g}{r^2}\right)-qg+g^2/2-\sigma \left(4{\left|g\right|}^2+2{\left|f\right|}^2\right)g=0\end{array}$ (19)

We used the solutions for S = 0 as our initial input to start our iterative process.

For a nonlinear problem like the one at hand, it has to keep in mind that stationary solutions only exist over a certain combinations and ranges of the system parameters. It is also possible that multiple solutions could exist.

We have obtained a set of stationary vortex solutions as shown in Figure 1, where |u|² is plotted over a {x,y}-plane. Figure 2 shows the contours of the same set of solutions. As a demonstration that stationary solutions have reached, Figure 3 shows that the pulse energy, E, for both u and v remain the same over the distance travelled, where

$\begin{array}{l}{E}_u={\displaystyle {\int }_{-20}^{20}{\displaystyle {\int }_{-20}^{20}{\left|u\right|}^2\text{d}x\text{d}y}}\\ E_v={\displaystyle {\int }_{-20}^{20}{\displaystyle {\int }_{-20}^{20}{\left|v\right|}^2\text{d}x\text{d}y}}\end{array}$ (20)

It should be noted that Equation (11) does not include a term for transmission loss. In reality, due to transmission loss both u and v would decrease in their intensities along z, needing periodic amplification to restore their pulse energy in long distance transmision.

5. Discussions

The stationary solutions for Equation (11) are complex electromagnetic waves. As they are complex quantities they exist only in abstract mathematics. In practice, we can only measure the absolute values of the pulse’s intensities or the pulse energies. But it is much easier to study electromagnetic waves, at least in this particular case, through the use of abstract complex algebra. No one would question whether complex electromagnetic waves do exist because they could not be measured directly. Therefore, a real event may be modelled using one or many abstract mathematical fields with abstract algebraic quantities as long they can be projected back to the real physical system that can be observed and measured.

It should be noted that Equation (11) contains terms associating with 1/r and this equation is used to cover a domain including the point at r = 0. These terms will become infinite unless the fields there have zero value. This is a necessary condition only due to the physical limitation that there is no infinite energy to support such infinite wave intensity. In a mathematical model, there is no such limitation as such a point could be replaced numerically by a nearby point where the field has a very large value. This limitation should be noted when a mathematical model is used to replace the physical event.

Equation (1) is much simpler than Equation (11). But both require integration along one variable, t, the time, in the case of Equation (1) and z, the distance travelled, in the case of Equation (11). We could solve Equation (1) analytically and there no problem to start integration forward from time zero or to start integration backward from a later time. The use of a different starting value could, however, give a completely different set of solutions. Equation (11) needs to be solved numerically. For a nonlinear problem like the one considered here, it is unlikely that using a solution for a later history could be used to recover the initial condition by integrating backward.

It should be noted that both equations are based on energy balance.

As with all other mathematical models, Equation (1) leads to Equation (3) in that two system parameters, c₁ and c₂ need to be found from observations or measurements, as we cannot find directly the material properties involved in Equation (1). Therefore, two sets of data are needed. Unless we can find the analytical solution, it is not possible to use one set of data and integrate forward or backward in time to discover the other set as reported in Reference [10].

For a more comprehensive model, more terms, or even equations, could be added to Equation (1). For example, an equation for Q could be developed to account for not just energy production but also energy changes due to expansion or gravity.

6. Conclusions

Our energy balance model has been used to show that the present observed background temperature of the universe need not be the remnant of a very high temperature past. The initial background temperature of the universe could be quite similar to what it is now.

Our model could accommodate an infinite initial temperature. This is done by following standard numerical technique: that is to replace the initial singular point by a very nearby point, for example, using time at 10⁻³⁵ second when temperature has assumed to have cooled to 10²⁷ K.

We have explained that an infinite temperature could exist in a mathematical model. In reality, this is not possible as such a temperature needs to be supported by an infinite amount of energy and this is physically not possible. Nevertheless, the universe’s huge mass, if converted, is near to infinite equivalent energy in mathematical terms; we are unable to explain, both mathematically and physically, from where this amount of energy could have come.