format: text

text:
# Adaptive sampling by information maximization

## Abstract

The investigation of input-output systems often requires a sophisticated choice of test inputs to make best use of limited experimental time. Here we present an iterative algorithm that continuously adjusts an ensemble of test inputs online, subject to the data already acquired about the system under study. The algorithm focuses the input ensemble by maximizing the mutual information between input and output. We apply the algorithm to simulated neurophysiological experiments and show that it serves to extract the ensemble of stimuli that a given neural system “expects” as a result of its natural history.

Biophysical systems often have many degrees of freedom and thus one needs large numbers of variables and parameters to describe them. Without strong prior knowledge about the intrinsic dynamics of such a system, one is left with inferring its function from data obtained by experiments or observations. Given a system where we control a set of “input” variables $`x=(x^{(1)},x^{(2)},\mathrm{},x^{(n)})`$ and measure another set of “output” variables $`y=(y^{(1)},y^{(2)},\mathrm{},y^{(m)})`$, we can actively manipulate the data acquisition by selecting the most informative test inputs. Yet how should one choose the test inputs to learn most about the input-output relation?

Within the classical Volterra-Wiener system identification methods Wiener (1958), the input space is sampled by drawing inputs from a probability distribution $`p(x)`$; a common choice is Gaussian “white noise”. However, not all aspects of the system’s input-output relation may be equally important. In neurobiology, for instance, one is especially interested in inputs $`x`$ about which a given sensory system conveys most information. In the spirit of importance sampling Itzykson and Drouffe (1989), one might therefore focus the data acquisition on those $`x`$ that contribute most to the information transfer. For a given input distribution, the information provided by a single input can then be quantified as $`I(x)=H_yH_y(x)`$ where $`H_y`$ is the entropy of the output distribution $`p(y)`$ and $`H_y(x)`$ is the entropy of the conditional probabilities $`p(y|x)`$ which characterize the input-output relation Shannon and Weaver (1949); Cover and Thomas (1991). Hence, the appropriate focusing is achieved by an input distribution $`p_{\text{opt}}(x)`$ that maximizes the mutual information $`I=I(x)`$ where the angular brackets denote averaging over $`p_{\text{opt}}(x)`$.

Without any information about the system and its input-output relation, the optimal input distribution $`p_{\text{opt}}(x)`$ is unknown. Any experimental test of the system must therefore start with drawing the test inputs from some predefined distribution $`p_\varphi (x)`$ that depends on a set of parameters $`\varphi =(\varphi ^{(1)},\mathrm{},\varphi ^{(L)})`$. Once data about the system has been acquired, however, one need not adhere to this initial choice of an input distribution. Instead, one should adapt the parameters or even the structure of $`p_\varphi (x)`$ to better focus on the important inputs. In this letter, we show how to systematically perform this adaptation. By iterating the adaptation procedure, the acquired data becomes ever more useful and the input distribution approaches the optimum.

Adapting the input distribution—For mathematical simplicity, we assume that both input and output take discrete values. Say that we have already tested the system with $`N`$ different inputs $`x_i`$ each of which was presented $`M_i`$ times while measuring the outputs $`y_{ij}`$ with $`i=1\mathrm{}N`$ and $`j=1\mathrm{}M_i`$. We define the set of all different output values measured so far by $`\{y_k:k=1\mathrm{}K\}`$. Our present knowledge about the system is summarized by the conditional probability that an output $`y_k`$ was obtained from the input $`x_i`$,

$$q(y_k|x_i)=\frac{1}{M_i}\underset{j=1}{\overset{M_i}{}}\delta _{y_{ij},y_k}.$$
(1)
The estimated probabilities $`q(y_k|x_i)`$ allow us to re-evaluate the relative importance of the inputs $`x_i`$ in terms of their potential contribution to the mutual information. To measure this contribution, we assign a probability or “weight” $`q(x_i)`$ to every input. Initially we assume that all inputs $`x_i`$ contribute equally and set $`q_1(x_i)=1/N`$. To find a combination of weights that maximizes the information transfer, we use the Blahut-Arimoto algorithm Arimoto (1972) and readjust the weights,

$$q_{n+1}(x_i)=\frac{1}{Z}q_n(x_i)\mathrm{exp}\left(\underset{k=1}{\overset{K}{}}q(y_k|x_i)\mathrm{log}\frac{q(y_k|x_i)}{q_n(y_k)}\right).$$
(2)
Here $`q_n(y_k)=_{i=1}^Nq(y_k|x_i)q_n(x_i)`$ and $`Z`$ is a normalization constant so that $`_{i=1}^Nq_{n+1}(x_i)=1`$. According to Eq. (2), the weight of an input $`x_i`$ is decreased if its conditional output distribution $`q(y_k|x_i)`$ is similar to the total output distribution $`q_n(y_k)`$. In contrast, the weight of an input $`x_i`$ is increased if the respective distributions differ. When Eq. (2) is iterated, the weights converge and reach a global maximum of the mutual information Arimoto (1972). In practice, we terminate the process once $`|1q_{n+1}(x_i)/q_n(x_i)|<ϵ`$ for all $`i`$ and some chosen precision $`ϵ`$ and set $`q_{\text{opt}}(x_i)=q_{n+1}(x_i)`$.

The weights or probabilities $`q_{\text{opt}}(x_i)`$ describe the relative frequencies with which the respective inputs $`x_i`$ should be drawn. Consequently, we need to adapt the parameters $`\varphi `$ so as to find a matching distribution $`p_\varphi (x)`$. Here we determine the new parameters $`\varphi `$ by maximizing the log-likelihood function Barlow (1989),

$$\mathrm{log}L(x_1,\mathrm{},x_N|\varphi )=\underset{i=1}{\overset{N}{}}q_{\text{opt}}(x_i)\mathrm{log}p_\varphi (x_i)$$
(3)
where the probabilities $`q_{\text{opt}}(x_i)`$ provide the appropriate weights. For some model distributions, e.g. Gaussians, the maximum can be found analytically. In general, however, one has to evaluate the maximum numerically.

The input distribution given by the new parameter values can be used to draw new test inputs, present them to the system and measure the respective outputs. After a certain amount of data has been acquired, the parameters $`\varphi `$ of the input distribution can be adapted again. The resulting iterative algorithm moves the input distribution towards an optimal ensemble.

Model quality and convergence—Every maximum of the mutual information with respect to $`p(x)`$ is a global maximum Cover and Thomas (1991). Hence, if the model distribution does not rule out any inputs, i.e., $`p_\varphi (x)>0`$ for all $`x`$ and $`\varphi `$, the estimates of the input-output relation, Eq. (1), converge, and therefore $`q_{\text{opt}}(x_i)p_{\text{opt}}(x_i)`$. Accordingly, the mutual information $`I_\text{D}=H_y^qH_y(x)_q`$ achieves the information capacity of the system; here, the index $`q`$ denotes that the respective quantities and averages are calculated with respect to $`q_{\text{opt}}(x_i)`$.

The model distribution $`p_\varphi (x)`$ converges towards an optimal fit of $`p_{\text{opt}}(x)`$. To control how well the model distribution captures the structure of the optimal distribution, one can check the mutual information achieved by the model, $`I_\text{M}=H_y^\varphi H_y(x)_\varphi `$, which is calculated with respect to $`r_\varphi (x_i)=p_\varphi (x_i)/[_{j=1}^Np_\varphi (x_j)]`$. The fraction $`\gamma `$ of the mutual information captured by the model is then defined as

$$\gamma =\frac{I_\text{M}}{I_\text{D}}$$
(4)
and provides a measure for the quality of the model. Hence, if $`\gamma `$ falls significantly below one, the model does no longer capture the structure of the optimal ensemble; in such a case, one might increase the complexity of the model.

In general, the algorithm will not be able to adapt the input ensemble if the presented inputs always result in the same output value. Similarly, there is no possibility to weight the inputs $`x_i`$ differently if every input elicits a new, different output. However, the latter problem can be solved by discretizing the output side into a smaller number of possible outputs. The input space, on the other hand, can be discretized as fine as needed without impeding the convergence of $`p_\varphi (x)`$.

Example—To illustrate the method, we study a numerical simulation of a Hodgkin-Huxley-type model neuron Wang and Buzsáki (1996). The model neuron transforms an input current $`I`$ into a voltage output $`V`$. For constant current values $`I<0\mu \text{A/cm}^2`$, the voltage approaches a stable equilibrium. For current values $`I>0\mu \text{A/cm}^2`$, the model undergoes a saddle-node bifurcation and generates periodically occurring action potentials, also called spikes Izhikevich (2000). Stochastic aspects of neural activity are incorporated by adding Gaussian white noise with a fixed standard deviation $`\sigma _\eta `$ and a cut-off frequency $`f_\eta `$ to the input.

We start with a simple one-dimensional parametrization of input and output. The inputs are 100-ms-long, discretized current steps ($`\mathrm{\Delta }I=1\mu \text{A/cm}^2`$), restricted to a physiologically realistic range of $`I=12\mathrm{}28\mu \text{A/cm}^2`$. The outputs are given by the number of spikes, $`C`$, during the corresponding time window. The resulting probabilistic relation of spike count versus current is displayed in Fig. 1(a).

For this one-dimensional input-output system, we can compute an exact solution of the information maximization problem. The optimal input distribution $`p_{\text{opt}}(I)`$ is depicted by the vertical bars in Fig. 1(b); the shape of $`p_{\text{opt}}(I)`$ corresponds to the slope of the input-output relation footnote (1). Note that there is a slight increase in the probabilities of inputs far below threshold ($`I10\mu \text{A/cm}^2`$). These inputs result almost certainly in a zero spike count output. At the same time, inputs closer to threshold ($`I9\mathrm{}1\mu \text{A/cm}^2`$) are more likely to produce spikes. As the optimal input distribution favors inputs that are more reliable, the inputs closer to threshold are neglected.

To study the performance of the iterative algorithm, we model the optimal input distribution by a truncated Gaussian. As initial parameter values, we choose a mean $`\varphi ^{(1)}=10\mu \text{A/cm}^2`$ and a standard deviation $`\varphi ^{(2)}=10\mu \text{A/cm}^2`$. In each iteration, we draw $`A`$ current values from the Gaussian, test them $`B`$ times on the system, and adapt the parameters. For our Gaussian model, the maximum likelihood estimate of the new parameters is given by $`\varphi ^{(1)}=_{i=1}^NI_iq_{\text{opt}}(I_i)`$ and $`\varphi ^{(2)}=\left[_{i=1}^N(I_i\varphi ^{(1)})^2q_{\text{opt}}(I_i)\right]^{1/2}`$.

The Gaussian model distributions are displayed in Fig. 1(b) for the first few iterations. Most of the current values drawn from the initial distribution fall below the spiking threshold of the neuron. Consequently, the algorithm shifts the Gaussian distribution into the spiking regime of the neuron. After about 10 iterations, the mutual information rate saturates at $`40`$ bits/sec. Since both the final Gaussian model $`p_\varphi (I)`$ and the optimal input distribution $`p_{\text{opt}}(I)`$ lead to approximately the same information transfer, the landscape of the mutual information with respect to the input distribution is relatively flat around the maximum; it suffices if the input distribution covers the relevant input range ($`I0\mathrm{}25\mu \text{A/cm}^2`$). Note, that due to the maximum-likelihood estimation, Eq. (3), the final Gaussian distribution has the same mean and variance as the optimal distribution.

Multi-dimensional example—The computational power of the algorithm becomes clearly visible for high-dimensional input spaces. As an example, consider the above model neuron when the input consists of time-varying, statistically stationary currents, discretized in time steps of $`\mathrm{\Delta }t_1`$. Following Strong et al. (1998), we slide overlapping windows of length $`T=n\mathrm{\Delta }t_1`$ across the input current trace and use the values within each window as input vector $`I_i=(I_i^{(1)},\mathrm{},I_i^{(n)})`$. For each of these inputs $`I_i`$, the output $`C_{ij}`$ is given by the spike times, discretized in time steps of $`\mathrm{\Delta }t_2=T/m`$, during the corresponding window. Hence, each input consists of $`n`$ real-valued numbers bounded within the interval $`I=12\mathrm{}28\mu \text{A/cm}^2`$, and each output consists of $`m`$ numbers whose values are either zero (no spike) or one (spike). Note, that we do not explicitly discretize the current values; we instead assume that every input $`I_i`$ is unique. For simplicity, we use a Gaussian input distribution. As the input is real and stationary, it suffices to use $`L=n/2+1`$ parameters for describing average and power spectrum of the current trace.

To test the system, we choose an initial distribution with an average $`\varphi ^{(1)}=0\mu \text{A/cm}^2`$ and a flat power spectrum with standard deviation $`\sigma =\left[_{i=1}^{n/2}\varphi ^{(i)}\right]^{1/2}=10\mu A/cm^2`$. For this prior, only 50% of the input values lie above threshold and the inputs will rarely lead to high firing rates, cf. Fig. 1. Consequently, we do not properly explore the full range of the input-output relation; if, for example, we test the system for 30 minutes with input currents drawn from this initial distribution, the information rate $`I_\text{D}`$ does not exceed $`300`$ bits/sec.

When using the iterative algorithm to adapt the parameters of the input ensemble, on the other hand, the information rate $`I_\text{D}`$ saturates around $``$ 670 bits/sec after about 20 minutes. Figure 2(a) shows how the power spectrum is shaped during the iterations. Only input frequencies below 500 Hz are well suited for the information transfer, the cut-off is roughly determined by the maximum firing rate of the model neuron. The overall increase in power leads to input currents that override the additive noise $`\eta `$ of the model neuron.

Initial conditions, convergence, and degeneracies—When the initial distribution is very narrow (flat power spectrum up to $`f_c=1000Hz`$, with $`\sigma =1\mu \text{A/cm}^2`$, $`\varphi ^{(1)}=20\mu A/cm^2`$, Fig. 2(b), dotted line), most of the input currents drive the neuron maximally and thereby very reliably. The strong initial bias leaves the algorithm with little maneuvering space for the parameter re-estimation so that it takes longer to approximate an optimal input distribution.

In the worst case, every input leads to the same output value. With the initial choice $`\varphi ^{(1)}=6\mu \text{A/cm}^2`$ and $`\sigma =10\mu \text{A/cm}^2`$, the input does not elicit any spikes during the first iterations, cf. Fig. 2(b), dashed line. However, once a spike has appeared, the statistics of the model distribution immediately moves into the direction of the statistics of the inputs $`I_i`$ that caused a spike. When the algorithm has tracked the relevant input range, a rapid increase of the information rate follows.

In the examples studied, the mutual information reaches approximately the same value independent of the initial conditions, cf. Fig. 2(b). Although there is always a clear preference for frequencies below 500 Hz, however, the parameters of the optimal input ensemble do not converge to the same set of values. Consequently, there is no unique combination of parameters that maximizes the mutual information; an observation that generalizes beyond the specific examples chosen. Nonetheless, the final input distributions always capture about $`\gamma =80`$% of the mutual information $`I_\text{D}`$, cf. Fig. 2(c).

In general, there might be “degenerate” subsets in stimulus space, i.e., sets of stimuli that lead to the same output value. In these cases, the total probability assigned to such a subset can be distributed in an arbitrary way on the subset and any statistical parameters $`\varphi `$ that depend on these subsets can assume different values without significant consequences for the information transfer.

Neurophysiological interpretation—Recent studies indicate that sensory neurons convey large amounts of information if the properties of the stimulus ensembles used match those of natural stimuli Rieke et al. (1995). Here we have shown how to extract a stimulus ensemble that conveys the maximum possible information without any prior knowledge. The proposed method could therefore serve to find the ensemble of stimuli that a given neuron naturally “expects”. Note that in contrast to previous online algorithms such as Alopex or Simplex Harth and Tzanakou (1974), we are not looking for a single optimal stimulus but rather for a complete ensemble of stimuli.

Our results demonstrate also that the optimal stimulus ensemble depends on the chosen criteria about what aspect of the output carries the relevant information. Hence, if the investigated model neuron conveys information in its average firing rate, then it best encodes slow-varying current values in the range $`I=0\mathrm{}23\mu \text{A/cm}^2`$. Synaptic inputs should therefore drive the neuron in the corresponding range. If, on the other hand, a neuron encodes its information in the precise timing of spikes, then the synaptic input should be of a more binary nature and either fully excite or fully inhibit the neuron. Measuring the time courses of a neuron’s membrane potential thus allows conclusions about the used neural code under optimal conditions.

###### Acknowledgements.

I thank A.V.M. Herz and M.B. Stemmler for stimulating discussions and H. Herzel for helpful comments on the manuscript. This work was supported by the DFG through the Innovationskolleg Theoretische Biologie and the Graduiertenkolleg 120.


text:
# Medium modification of the nucleon axial form factor

## Abstract

We study the modification of the nucleon axial form factor in nuclear matter. The internal quark substructure of the nucleon is self-consistently described by the quark meson coupling model. We find that the axial form factor of the bound nucleon is quenched considerably from that of the free nucleon. The axial vector coupling constant, $`g_A`$, is reduced by roughly 10% at normal nuclear matter density and the axial form factor varies within 8% for moderate momentum transfer.

ADP-01/52-T484

There is strong evidence that hadron properties must undergo substantial modifications in nuclear medium . A number of experiments, such as the variation of nucleon structure functions in lepton deep-inelastic scattering off nuclei (the nuclear EMC effect) , the quenching of the axial vector coupling constant, $`g_A`$, in nuclear $`\beta `$-decay , and the missing strength of the response functions in nuclear quasielastic electron scattering , have stimulated investigations of whether or not quark degrees of freedom play any significant role. Though the conventional interpretation arising through polarization effects and other hadronic degrees of freedom ($`\mathrm{\Delta }`$-excitations, meson exchange currents, etc.) undoubtedly play a role , it is rather interesting to explore, as well, the possible effects of a change in the internal structure of the bound nucleon.

The successes of quantum hadrodynamics (QHD) leave little doubt that relativistic nuclear phenomenology is essential in describing the bulk properties of nuclear matter as well as the properties of finite nuclei. To incorporate the substructure of the nucleon in a relativistic nuclear framework, Guichon proposed a successful hybrid model (the quark-meson coupling (QMC) model) , where nuclear matter is described in terms of non-overlapping, MIT bag nucleons . The model was later developed to describe finite nuclei as well as the properties of other hadrons in medium . By analogy with QHD, QMC describes the bulk properties of nuclear systems using scalar $`(\sigma )`$ and vector $`(\omega )`$ meson mean fields. However, the nucleon bound in a nuclear medium here is no longer a point-like particle, it has substructure – quarks confined inside the nucleon bag. It is the quark, rather than the nucleon itself, that is coupled to the $`\sigma `$ and $`\omega `$ fields directly. As a result, the internal structure of the bound nucleon is modified by the surrounding medium with respect to that in free space .

Within QMC, the small mass of the quark implies that the lower component of the quark wave function will be enhanced rapidly §§§Considerably more rapidly than that of the nucleon in QHD and for smaller values of the mean-fields. by the change of its environment (as the $`\sigma `$ field strength increases), with a consequent decrease in the scalar baryon density. As the scalar baryon density itself is the source of the $`\sigma `$ field, this provides a mechanism for the saturation of nuclear matter where the quark substructure plays a vital role. The extra degrees of freedom, corresponding to the internal structure of the nucleon, lead to a reasonable value for the nuclear incompressibilty, once the corresponding quark and meson coupling constants, $`g_\sigma ^q`$ and $`g_\omega ^q`$, are determined to reproduce the empirical values for the saturation density and binding energy of symmetric nuclear matter.

In our previous work, we predicted medium modifications of nucleon electromagnetic form factors . Such a medium effect appears to be supported by a recent experiment at Mainz , which measured the polarization transfer in the $`{}_{}{}^{4}He(\stackrel{}{e},e^{}\stackrel{}{p})`$ reaction. The polarization transfer double ratio, $`(P_x^{}/P_z^{})_{\mathrm{He}}/(P_x^{}/P_z^{})_{\mathrm{free}}`$, tends to favor the RDWIA calculations using a medium modified proton form factor. Further study of this reaction has been carried out at Jefferson Lab (E93-049)  and other related experiments have been proposed. In this paper, we study the axial form factor of the bound nucleon in symmetric nuclear matter. This is of particular interest in the light of plans to build very high intensity neutrino beams in the near future.

The Lagrangian density for the QMC model is

$`_q`$ $`=`$ $`\overline{q}(i\gamma ^\mu _\mu m_q)q\theta _VB_0\theta _V+g_\sigma ^q\overline{q}q\sigma g_\omega ^q\overline{q}\gamma _\mu q\omega ^\mu {\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+{\displaystyle \frac{1}{2}}m_\omega ^2\omega ^2,`$ (1)
where $`m_q`$ is the current quark mass, $`B_0`$ is the bag constant in vacuum, $`g_\sigma ^q`$ and $`g_\omega ^q`$ are the corresponding quark and meson coupling constants, and $`\theta _V`$ is a step function which is one inside the bag volume and zero outside.

In mean-field-approximation, the meson fields are treated as classical fields, and the quark field $`q(x)`$ inside the bag satisfies the equation of motion

$$[i\gamma ^\mu _\mu (m_qg_\sigma ^q\overline{\sigma })g_\omega ^q\overline{\omega }\gamma ^0]q(x)=0,$$
(2)
where $`\overline{\sigma }`$ and $`\overline{\omega }`$ denote the constant mean values of the scalar and the time component of the vector field in symmetric nuclear matter. The normalized solution for the lowest state of the quark is given by

$$q(t,\stackrel{}{r})=\frac{N_0}{\sqrt{4\pi }}e^{iϵ_qt/R}\left(\begin{array}{c}g(r)\\ i𝝈\widehat{𝐫}f(r)\end{array}\right)\theta (Rr)\chi _q,$$
(3)
where

$`g(r)`$ $`=`$ $`j_0(xr/R),f(r)=\beta _qj_1(xr/R),`$ (4)
$`ϵ_q`$ $`=`$ $`\mathrm{\Omega }_q+g_\omega ^q\overline{\omega }R,\beta _q=\sqrt{{\displaystyle \frac{\mathrm{\Omega }_qm_q^{}R}{\mathrm{\Omega }_q+m_q^{}R}}},`$ (5)
$`N_0^2`$ $`=`$ $`2R^3j_0^2(x)[\mathrm{\Omega }_q(\mathrm{\Omega }_q1)+m_q^{}R/2]/x^2,`$ (6)
with $`\mathrm{\Omega }_q\sqrt{x^2+(m_q^{}R)^2}`$, $`m_q^{}m_qg_\sigma ^q\overline{\sigma }`$, $`R`$ ($`R_0`$) the bag radius (free space), and $`\chi _q`$ the quark Pauli spinor. The eigen-frequency, $`x(x_0)`$, of this lowest mode in medium (free space) is determined by the boundary condition at the bag surface,

$$j_0(x)=\beta _qj_1(x).$$
(7)
The form of the quark wave function in Eq. (3) is almost identical to that of the solution in free space. However the parameters in the expression have to be substantially modified by the surrounding nuclear medium. Note that as the value of $`g_\sigma ^q\overline{\sigma }`$ is usually much larger than $`m_q`$, the quantity, $`\beta _q`$, becomes larger than unity, which means that the lower component of the Dirac spinor is enhanced. In other words, the quarks in the nucleon embedded in the nuclear medium are more relativistic than those in a free nucleon.

The mean values of the scalar ($`\overline{\sigma }`$) and vector ($`\overline{\omega }`$) fields in symmetric nuclear matter are self-consistently determined by solving the following set of equations:

$`\overline{\omega }`$ $`=`$ $`{\displaystyle \frac{g_\omega \rho }{m_\omega ^2}},`$ (8)
$`\overline{\sigma }`$ $`=`$ $`{\displaystyle \frac{g_\sigma }{m_\sigma ^2}}C(\overline{\sigma })\rho _s={\displaystyle \frac{g_\sigma }{m_\sigma ^2}}C(\overline{\sigma }){\displaystyle \frac{4}{(2\pi )^3}}{\displaystyle ^{k_F}}d^3k{\displaystyle \frac{m_N^{}(\overline{\sigma })}{\sqrt{m_N^2(\overline{\sigma })+k^2}}},`$ (9)
$`m_N^{}(\overline{\sigma })`$ $`=`$ $`{\displaystyle \frac{3\mathrm{\Omega }_q(\overline{\sigma })}{R}}{\displaystyle \frac{z_0}{R}}+{\displaystyle \frac{4}{3}}\pi R^3B_0,`$ (10)
$`{\displaystyle \frac{m_N^{}(\overline{\sigma })}{R}}`$ $`=`$ $`0,`$ (11)
where $`\rho `$ ($`\rho _s`$) is the baryon (scalar) density and $`k_F`$ is the nucleon Fermi momentum, $`g_\sigma =3g_\sigma ^qS(0)`$, $`g_\omega =3g_\omega ^q`$, and the quantity, $`C(\overline{\sigma })`$, is defined by

$$C(\overline{\sigma })S(\overline{\sigma })/S(0)=\left(\frac{m_N^{}(\overline{\sigma })}{\overline{\sigma }}\right)/g_\sigma ,$$
(12)
with $`S(\overline{\sigma })=_{\mathrm{bag}}d^3r\overline{q}(\stackrel{}{r})q(\stackrel{}{r})`$. Using the quark wave function for the MIT bag, Eq.(3), $`S(\overline{\sigma })`$ can be explicitly evaluated:

$$S(\overline{\sigma })=[\mathrm{\Omega }_q/2+m_q^{}R(\mathrm{\Omega }_q1)]/[\mathrm{\Omega }_q(\mathrm{\Omega }_q1)+m_q^{}R/2].$$
(13)
The second term in Eq.(10), ($`z_0/R`$), has multiple roles and it parametrizes the sum of the zero point energy, gluon corrections, and the part of the center-of-mass (c.m.) motion.

The process to solve this coupled system is as follows: we first determine $`z_0`$ and $`B_0`$ by requiring the free nucleon mass to be $`m_N(\rho =0)=939`$ MeV and by imposing the stability condition, Eq.(11), for a given bag radius, $`R_0`$ (treated as an input parameter). After that, we solve the coupled set of equations at normal nuclear matter density, $`\rho _0`$, and determine the coupling constants, $`g_\sigma ^q`$ and $`g_\omega ^q`$, required to reproduce nuclear saturation properties. With these parameters, we can solve the whole problem for each finite nuclear matter density, $`\rho `$, self-consistently. Typically the quark r.m.s. radius, $`r_q^{}`$, calculated by the bag wave function is slightly increased, although the bag radius, $`R`$, decreases by a few percent at normal nuclear matter density. The properties of these self-consistent solutions as a function of the nuclear matter density, a typical set of parameters and the value of the nuclear incompressibility, $`K`$, as well as the possible medium dependence of the bag constant, can be seen in Ref. .

It is worthwhile to note that the self-consistency condition in QMC is identical to that in QHD, except that in QHD one has $`C(\overline{\sigma })=1`$ in Eq.(9, which corresponds to a point-like nucleon. Within QMC all information on the internal structure of the nucleon is contained in $`C(\overline{\sigma })`$.

Once the quark wave function in the bound nucleon is determined, one can proceed to calculate the nucleon axial form factors, which are defined as follows:

$$p^{}s^{}\left|A_a^\mu (0)\right|ps=\overline{u}_s^{}(p^{})\left[G_A(Q^2)\gamma ^\mu +\frac{G_P(Q^2)}{2m_N}(p^{}p)^\mu \right]\gamma _5\frac{\tau _a}{2}u_s(p),$$
(14)
where $`Q^2(p^{}p)^2`$, $`u_s(p)`$ is the nucleon Dirac spinor. Here we shall focus on the axial vector form factor, $`G_A(Q^2)`$, since the induced pseudoscalar form factor, $`G_P(Q^2)`$, is dominated by the pion pole and thus can be derived using the familiar PCAC relation . The relevant axial current operator is then simply

$`A_a^\mu (x)`$ $`=`$ $`{\displaystyle \underset{f}{}}\overline{q}_f(x)\gamma ^\mu \gamma _5{\displaystyle \frac{\tau _a}{2}}q_f(x)\theta (Rr),`$ (15)
where $`q_f(x)`$ is the quark field operator for flavor $`f`$.

In order to remove the spurious c.m. motion, we construct the momentum eigenstate of a baryon via the Peierls-Thouless (PT) projection method,

$$\mathrm{\Psi }_{\mathrm{PT}}(\stackrel{}{x}_1,\stackrel{}{x}_2,\stackrel{}{x}_3;\stackrel{}{p})=N_{\mathrm{PT}}e^{i\stackrel{}{p}\stackrel{}{x}_{\mathrm{c}.\mathrm{m}.}}q(\stackrel{}{x}_1\stackrel{}{x}_{\mathrm{c}.\mathrm{m}.})q(\stackrel{}{x}_2\stackrel{}{x}_{\mathrm{c}.\mathrm{m}.})q(\stackrel{}{x}_3\stackrel{}{x}_{\mathrm{c}.\mathrm{m}.}),$$
(16)
where $`N_{\mathrm{PT}}`$ is a normalization constant, $`\stackrel{}{p}`$ the total momentum of the baryon, and $`\stackrel{}{x}_{\mathrm{c}.\mathrm{m}.}=(\stackrel{}{x}_1+\stackrel{}{x}_2+\stackrel{}{x}_3)/3`$ is the center of mass of the baryon (we assume equal mass quarks here). It can be shown that the PT wave function satisfies the condition of translational invariance. Using Eqs. (15) and (16), the nucleon axial form factors can be expressed as

$`G_A(Q^2)`$ $`=`$ $`{\displaystyle \frac{5}{3}}{\displaystyle d^3r\left\{\left[g^2(r)f^2(r)\right]j_0(Qr)+2f^2(r)\frac{j_1(Qr)}{Qr}\right\}K(r)/D_{\mathrm{PT}}},`$ (17)
$`D_{\mathrm{PT}}`$ $`=`$ $`{\displaystyle d^3r\rho _q(r)K(r)},`$ (18)
where $`D_{\mathrm{PT}}`$ is the normalization factor, $`\rho _q(r)g^2(r)+f^2(r)`$, and $`K(r)=d^3z\rho _q(\stackrel{}{z})\rho _q(\stackrel{}{z}\stackrel{}{r})`$ is the recoil function which accounts for the correlation of the two spectator quarks.

At this stage, there is no satisfactory covariant treatment for the MIT bag model. On the other hand, relativistic effects are important for most dynamic variables, especially for form factors at large momentum transfer. They lead to a sizable correction for the r.m.s. radius of the nucleon. In this paper, we use a semi-phenomenological method to account for the relativistic corrections which is consistent with a mean field treatment of the QMC model. Since a static MIT bag is an extended spherical object, it would be deformed if it were viewed in a moving frame of reference. It is crucial to include this Lorentz contraction of the bag for calculating form factors at moderate momentum transfer. In the prefered Breit frame, the resulting form factors can be expressed through a simple rescaling, i.e.,

$`G_A(Q^2)`$ $`=`$ $`\left({\displaystyle \frac{m_N^{}}{E^{}}}\right)^2G_A^{\mathrm{sph}}(Q^2m_{N}^{}{}_{}{}^{2}/E_{}^{}{}_{}{}^{2}),`$ (19)
where $`E^{}=\sqrt{m_{N}^{}{}_{}{}^{2}+Q^2/4}`$ and $`G_A^{\mathrm{sph}}(Q^2)`$ are the form factors calculated with the spherical bag wave function. The scaling factor in the argument arises from the coordinate transformation of the struck quark whereas the prefactor, $`(m_N^{}/E^{})^2`$, comes from the reduction of the integral measure of two spectator quarks in the Breit frame. The axial radius squared is given by,

$$r_A^2=\frac{6}{g_A}\frac{dG_A(Q^2)}{dQ^2}|_{Q^20},$$
(20)
which results in $`12/m_A^2`$ for a dipole form: $`G_A(Q^2)=1/(1+Q^2/m_A^2)^2`$.

Note that the pion cloud of the nucleon only plays an indirect role in calculating the axial vector form factor, in contrast to the electromagnetic properties of the nucleon. In fact, the matrix element of the pionic axial current, $`f_\pi ^\mu \pi `$, vanishes in any chiral bag model if the pion cloud exists in all space .

The nucleon axial form factor in free space is illustrated in Fig. 1. The experimental data from neutrino scattering is not shown directly, rather we show a dipole form with the mass parameter in the range found by Kitagachi et al. . (We note that the pion electroproduction data  is consistent with this after the small correction from chiral perturbation theory is applied .) The corrections for the center-of-mass motion and Lorentz contraction lead to a significant improvement over static bag model calculations, in particular, at moderate momentum transfers . The results are not sensitive to the current quark mass, so we use $`m_q=0`$ in this paper. The axial radii are 0.587, 0.614 and 0.640 fm, corresponding to the bag radii of 0.90, 0.95 and 1.00 fm, respectively. These r.m.s. radii are in good agreement with the experimental value $`(0.635\pm 0.023)`$ fm . The axial coupling constant, $`g_AG_A(Q^2=0)`$, is about 1.14 after the c.m. correction, which is about 5% increase from the static MIT bag value, 1.09. Finite quark mass and pion renormalization may be expected to lead to further corrections at the level of 10%.

Fig. 2 shows the medium dependence of the axial form factors at several momentum transfers (with R0 = 0.95 fm). For $`Q^2<2`$ GeV2, the axial form factor decreases as the density $`\rho `$ increases. The axial vector coupling constant (the solid line in this figure), $`g_A`$, would be reduced by 12% at normal nuclear matter density, $`\rho _0`$, and 9% at the average density of finite nuclei, 0.7$`\rho _0`$, where the latter is to be compared with the experimental 20% reduction seen in Gamow-Teller transitions .

The momentum dependence of the axial form factors (R0 = 0.95 fm again) is shown in Fig. 3. The nuclear medium does indeed modify the shape of the momentum dependence. The medium effect increases as the density increases and tends to be larger for smaller momentum transfers. However the overall effect is less than roughly 8%, depending on the bag radius.

In conventional nuclear physics, the nucleon is immutable and the electroweak properties of a nucleus are often described by a combination of individual nucleon contributions and various meson exchange current corrections. This is particularly important for the axial charge of finite nuclei . For the quenching of $`g_A`$, conventional mechanisms (medium polarization) also exist, however it is fundamentally interesting to explore this new mechanism arising from possible changes of the internal structure of the nucleon. Clearly it would be very valuable to study the effect of the medium modified axial form factor in the context of neutrino-nucleus scattering and the solar neutrino problem .

In summary, the nucleon axial form factor is significantly modified in nuclear matter through a new mechanism, namely the change of the internal substructure of the bound nucleon. We have calculated the density dependence of the axial form factor of the bound nucleon in the QMC model and found that it is quenched considerably compared with that of the free nucleon. The axial vector coupling constant, $`g_A`$, is reduced by roughly 10% at $`\rho _0`$ and the axial form factor would vary within 8% for moderate momentum transfers.

D.H.Lu is grateful to the Y.C.Tang Disciplinary Development Fund in Zhejiang University and would like to acknowledge the warm hospitality of the CSSM in Adelaide University. This work was supported in part by the National Natural Science Fund of China and by the Australian Research Council.


text:
# Trans-Planckian Physics and the Spectrum of Fluctuations in a Bouncing Universe

## I Introduction

The dependence of the spectrum of cosmological fluctuations in an inflationary Universe on hidden assumptions about the physics on length scales much smaller than the Planck length has been recently studied in Refs. BM1 ; MB2 ; N1 . A priori this dependence comes about since in typical scalar-field-driven inflationary models the duration of the period of expansion is so long that the physical wavelengths of comoving modes which correspond to the present large-scale structure of the Universe are much smaller than the Planck length at the beginning of inflation. In weakly coupled scalar field models of inflation, the spectrum of fluctuations is calculated by assuming that fluctuation modes start out in the vacuum (i.e. minimum energy density) state at the beginning of inflation and subsequently evolve as determined by the equations of motion. It is hence not unreasonable to expect that modifications of the physics on trans-Planckian scales could affect the final spectrum of fluctuations.

While the correct theory of trans-Planckian physics is not known, possible effects of the new physics can be modeled by changes in the dispersion relation of the fields corresponding to linear cosmological fluctuations. Such an approach was used earlier Unruh ; CJ  to study the possible dependence of the spectrum of black hole radiation on trans-Planckian physics. Indeed, for a class of dispersion relations which deviate so strongly from the usual one such that the evolution of the mode functions is highly non-adiabatic and the effective frequency becomes imaginary, it was found that the spectrum of fluctuations is modified in a significant way compared to what is obtained for a linear dispersion relation BM1 ; MB2 . As shown in MB3 ; NP2 ; Starob , the final spectrum of fluctuations is unchanged if the change in the dispersion relation is not too drastic, the wave function evolves adiabatically, and hence the WKB approximation is valid throughout the evolution.

Two major deficiencies in the previous work concerned the choice of initial conditions and the fact that complex dispersion relations were required in order to obtain significant deviations from the standard results. In the previous papers, the assumption was made that the state starts out as the one minimizing the energy density. However, if the dispersion relation differs greatly from the linear one, in particular if it becomes complex, the physical motivation for choosing this state becomes unclear. This problem does not arise in a bouncing Universe. In this case, modes corresponding to present cosmological scales could well have had a physical wavelength smaller than the Planck length at the bounce point, but at early (pre-bounce) times, the wavelength was larger than the Planck length. In this context, initial conditions can be set up at some very early pre-bounce time during which the mode obeys a linear dispersion relation and the physical meaning of the minimum energy density state is clear. It is thus interesting to study the spectrum of fluctuations in a bouncing Universe which is asymptotically flat such that the choice of initial vacuum becomes unambiguous. It is also interesting to exhibit a case where complex dispersion relations are not needed. Although complex frequencies are standard in classical physics and in quantum mechanics, they seem to be more problematic in the context of quantum field theory.

In this article, we find an example of a bouncing Universe in which the spectrum of fluctuations depends sensitively on changes in the physics at length scales smaller than the Planck length, without requiring complex frequencies at high wavenumbers, thus improving on the second major deficiency of previous work. Another reason for investigating the dependence of the spectrum of fluctuations on modifications of trans-Planckian physics in a bouncing Universe is that one could expect (however, as we show, this is not correct) that the deviations in the spectrum incurred during the pre-bounce period when the wavelength is blue-shifting in the trans-Planckian regime will cancel with the deviation in the post-bounce trans-Planckian regime when the mode is being red-shifted. Such a result would then lead to the expectation that the results on the dependence of the spectrum of fluctuations on changes in the trans-Planckian physics might be different in the case of black holes and (non-bounce) inflationary cosmology. A third reason for analyzing a bouncing Universe is that there has recently been a lot of interest in the field of string cosmology in bouncing Universe models. In the Einstein frame, the pre-big-bang scenario PBB  starts in a collapsing dilaton-dominated phase, and the same is true in the Ekpyrotic scenario KOST . Our arguments show that trans-Planckian effects could change the predictions of standard cosmological perturbation theory in these examples 333Note that there is disagreement on the result of the linear theory of fluctuations in the Ekpyrotic scenario KOST ; Lyth ; BF ; KOST2 ; Hwang ; MPPS .

A model for the evolution of the scale factor in a bouncing inflationary Universe is $`a(t)=a_0\mathrm{cosh}(Ht)`$, where $`a_0`$ denotes the value of $`a(t)`$ at the center of the bounce, and $`H`$ is the Hubble expansion (contraction) rate during the period of exponential expansion (contraction). The idea of the calculation is to assume an initial spectrum of fluctuations at some early time $`t_\mathrm{f}`$ when the wavelength of all modes of interest is much larger than the Planck length, to propagate the perturbation modes through the bounce and calculate the final spectrum at time $`t_\mathrm{f}`$. In this model, one of the difficulties encountered in the previous work BM1 ; MB2 , namely the problem of specifying the initial state in the trans-Planckian regime, is overcome. However, it is necessary now to set up initial conditions on wavelengths which typically also are much larger than the Hubble radius. The latter problem can be overcome by considering a bouncing Universe which is radiation-dominated in the asymptotic past and future, e.g. AS ; BD  $`a(\eta )=\sqrt{\alpha ^2+\beta ^2\eta ^2}`$, where $`\alpha `$ and $`\beta `$ are constants. Another toy model in which it is possible to set up well motivated initial conditions on scales much smaller than the Hubble radius and at the same time much larger than the Planck scale is a model in which the Universe is asymptotically flat both at large positive and negative times:

$$a(\eta )=\mathrm{}_0\frac{\mathrm{}_0\mathrm{}_\mathrm{b}}{1+(\eta /\eta _0)^2}$$
(1)
where $`\mathrm{}_0`$ and $`\mathrm{}_\mathrm{b}`$ (with $`\mathrm{}_\mathrm{b}\mathrm{}_0`$) denote the asymptotic size of the Universe and the size at the bounce, respectively, and $`\eta _0`$ determines the time scale of the bounce, i.e. the time over which the scale factor changes. In the above, we are using conformal time $`\eta `$, and employing the convention that the scale factor carries dimension of length while the comoving coordinates (including conformal time) are dimensionless.

Finally, a comment is in order on the issue of backreaction. It has been shown in Refs. Tanaka ; Starob  that the energy density of the particles created by trans-Planckian physics can contribute in a significant way to the background energy density. We do not address this interesting question since this clearly comes as a second step once a case where the spectrum is modified has been explicitly exhibited. In this paper we concentrate on the first question since it would be useless to treat the backreaction question if no well-motivated example where the spectrum is modified can be found. However, it should be clear that, once this is done (as shown in the present article), this problem becomes the central question in the study of trans-Planckian physics in cosmology.

This article is organized as follows. In section II, we review the arguments which show that no change in the final spectrum of fluctuations of a scalar matter field in a bouncing Universe is expected independent of the dispersion relation, provided that the WKB approximation remains justified. We give a qualitative reason to expect changes in the amplitude of the final spectrum for dispersion relations for which the adiabatic approximation breaks down. As mentioned above, our analysis assumes that perturbation theory remains justified and back-reaction can be neglected Tanaka . In Section III we study a concrete model, namely the asymptotically flat bouncing Universe given by (1) with a dispersion relation modified in the trans-Planckian regime according to the prescription of Corley and Jacobson CJ . We show that, as expected, both the overall amplitude and the shape of the spectrum differ from what is obtained for the linear dispersion relation, without ever requiring the effective frequency to become imaginary. The change in the spectrum is produced by an interesting combination of effects due to the modified dispersion relation and the driving term for Parker particle production Parker .

## II Method and Qualitative Considerations

The equation of motion for a minimally coupled free scalar matter field $`\mathrm{\Phi }(\eta ,𝒙)`$ in a Universe with scale factor $`a(\eta )`$ in momentum space takes on the simple form

$$\mu ^{\prime \prime }+\left(n^2\frac{a^{\prime \prime }}{a}\right)\mu =\mathrm{\hspace{0.17em}0},$$
(2)
making use of the re-scaling

$$\mathrm{\Phi }(\eta ,𝒙)=\frac{1}{(2\pi )^{3/2}}\frac{1}{a(\eta )}d𝒏\mu (n,\eta )e^{i𝒏𝒙},$$
(3)
where $`𝒏`$ denotes the comoving wavenumber linked to the physical wavenumber by the relation $`𝒌=𝒏/a(\eta )`$. In the above, a prime denotes the derivative with respect to conformal time. The previous equations are valid for a spacetime with flat spatial sections. The advantage of considering a test scalar field and/or gravitational waves is that we do not need to address the origin of the dynamics of the scale factor. This is in contrast to the case of density perturbations. However, the Friedmann equation reads $`^2/a^2+K/a^2=\kappa _i\rho _i/3`$, where $`a^{}/a`$ and $`\rho _i>0`$ is the energy density. This means that a bounce is consistent with the Einstein equations only if $`K=1`$. For $`K=0`$, it would imply $`\rho =0`$ for each component. Therefore, it seems that it is inconsistent to use Eq. (2) if the scale factor behaves such that there is a bounce. There are two ways out of the previous argument. The first one is the following. If $`K=1`$, then Eq. (2) takes the form $`\mu ^{\prime \prime }+(n^2Ka^{\prime \prime }/a)\mu =0`$. In addition the eigenfunctions of the three-dimensional Laplacian operator are no longer planes waves. However, if we restrict ourselves to modes with wavelengths much smaller than the curvature radius, then they do not feel the curvature of the spacelike section and we can safely work with Eq. (2). It could be checked that this is indeed the case for the wavenumbers considered in the example presented in the following section. Another possibility is to use the “nonsingular Universe construction” bms . In this case, it is definitely possible to get a spatially flat bounce. What happens is that the higher derivative terms become important at the bounce and enable the evasion of the previous argument by supplying other terms in the analog of the first Friedmann equation. In conclusion, although in this article we consider a scale factor with a bounce, we can nevertheless work with Eq. (2), i.e. we let the nonsingular Universe construction enter only via the scale factor.

The method introduced in Unruh ; CJ  to study the dependence of the spectrum of fluctuations on trans-Planckian physics is to replace the linear dispersion relation $`\omega _{_{\mathrm{𝐩𝐡𝐲𝐬}}}=k`$ by a non standard dispersion relation $`\omega _{_{\mathrm{𝐩𝐡𝐲𝐬}}}=\omega _{_{\mathrm{𝐩𝐡𝐲𝐬}}}(k)`$ where the function $`\omega _{_{\mathrm{𝐩𝐡𝐲𝐬}}}(k)`$ is a priori arbitrary. In the context of cosmology, it has been shown in Ref. BM1 ; MB2  that this amounts to replacing $`n^2`$ appearing in (2) with $`n_{\mathrm{eff}}^2(n,\eta )`$ defined by

$$n^2n_{\mathrm{eff}}^2(n,\eta )a^2(\eta )\omega _{_{\mathrm{𝐩𝐡𝐲𝐬}}}^2\left[\frac{n}{a(\eta )}\right].$$
(4)
For a fixed comoving mode, this implies that the dispersion relation becomes time-dependent. Various forms for $`n_{\mathrm{eff}}`$ have been used, e.g.

$$n_{\mathrm{eff}}(n,\eta )=n\frac{\lambda (\eta )}{\mathrm{}__\mathrm{C}}\mathrm{tanh}^{1/p}\left[\left(\frac{\mathrm{}__\mathrm{C}}{\lambda (\eta )}\right)^p\right],n_{\mathrm{eff}}^2(n,\eta )=n^2+n^2b_m\left[\frac{\mathrm{}__\mathrm{C}}{\lambda (\eta )}\right]^{2m},$$
(5)
the first being a generalization of the relation studied by Unruh Unruh  whereas the second one is a generalization of the relations used by Corley and Jacobson CJ . In the above, $`\lambda (\eta )`$ is the physical wavelength, $`\mathrm{}__\mathrm{C}`$ is the cutoff length (which can be taken to be the Planck length), $`p`$ and $`m`$ are integers, and $`b_m`$ is a real number which can be either positive or negative. Thus, the equation of motion to be analyzed is

$$\mu ^{\prime \prime }+\left[n_{\mathrm{eff}}^2(n,\eta )\frac{a^{\prime \prime }}{a}\right]\mu =\mathrm{\hspace{0.17em}0}.$$
(6)
For fixed comoving wavenumber $`n`$, the evolution of $`\mu `$ depends crucially on which region $`\eta `$ lies in. The first region is the trans-Planckian region with

$$\lambda (\eta )\mathrm{}__\mathrm{C}(\mathrm{Region}1).$$
(7)
Let us first consider the standard inflationary scale factor and a monotonic dispersion relation \[like the first one in Eq. (5) or the second one with $`b_m>0`$\]. This means that, initially, the term $`n_{\mathrm{eff}}^2(n,\eta )`$ dominates in Eq. (6). The initial conditions are fixed in this region and since the WKB approximation is applicable, we can choose the initial state as the “minimizing energy state” Brown78 . Then, the (positive frequency) solution is given by

$$\mu (\eta )\frac{1}{\sqrt{2n_{\mathrm{eff}}(n,\eta )}}\mathrm{exp}\left[i_{\eta _\mathrm{i}}^\eta n_{\mathrm{eff}}d\tau \right],$$
(8)
where $`\eta _i`$ is some initial time.

The second region corresponds to

$$\lambda \mathrm{}__\mathrm{C}\mathrm{and}n^2a^{\prime \prime }/a(\mathrm{Region}2).$$
(9)
In this region the mode has reached the linear part of the dispersion relation. The general solution in Region 2 is the plane wave

$$\mu =B_1\mathrm{exp}(in\eta )+B_2\mathrm{exp}(in\eta )$$
(10)
with constant coefficients $`B_1`$ and $`B_2`$. For the standard dispersion relation, the initial conditions are fixed in this region. The usual choice of the vacuum state is $`B_1=1/\sqrt{2n}`$, $`B_2=0`$. In general $`B_1`$ and $`B_2`$ are determined by the matching conditions between Regions 1 and 2. However, if the dynamics is adiabatic throughout (in particular if the $`a^{\prime \prime }/a`$ term is negligible), the WKB approximation holds and the solution is always given by (8). Therefore, if we start with $`B_2=0`$ and uses this solution, one finds that $`B_2`$ remains zero at all times. Deep in the region where $`n_{\mathrm{eff}}n`$ the solution becomes

$$\mu (\eta )\frac{1}{\sqrt{2n}}\mathrm{exp}(i\varphi in\eta ),$$
(11)
i.e. the standard vacuum solution times a phase which will disappear when we calculate the modulus. The phase $`\varphi `$ is given by $`\varphi _{\eta _\mathrm{i}}^{\eta _1}n_{\mathrm{eff}}d\tau `$, where $`\eta _1`$ is the time at which $`n_{\mathrm{eff}}n`$.

The situation changes dramatically if we consider non-monotonic dispersion relations. This is the case if $`b_m<0`$ in Eq. (5) or for the dispersion relation introduced in Ref. Mersini . Two new features can occur. Either the dispersion relation can become complex (this is in general the case if $`b_m<0`$) or the term $`a^{\prime \prime }/a`$ dominates (this is the case for the dispersion relation introduced in Ref. Mersini  or for the Corley/Jacobson dispersion relation with $`b_m<0`$ for a certain range of comoving wavenumbers). In this case, the WKB approximation is violated in Region 1 and we expect changes in the final spectrum. Unfortunately, in the context of the standard increasing inflationary scale factor, one also looses the ability to fix natural initial conditions. However, this is no longer true if the spacetime is asymptotically flat because then, at infinity, the term $`a^{\prime \prime }/a`$ goes to zero, see also Ref. LLMU . Therefore, in this case, we can choose well-motivated initial conditions. A good example of such a situation is provided by a bouncing Universe. From the above considerations, we expect that in this case the final spectrum is modified and the initial conditions can be fixed naturally. As an additional benefit, complex frequencies can be avoided. In the next section, we consider a toy model where these qualitative arguments can be implemented concretely, at the level of equations.

## III A Specific Example

We will now illustrate the qualitative arguments of the previous section with a concrete quantitative example. We take the asymptotically flat bouncing Universe given by (1). We consider the type of modified dispersion relation for which in the case of an expanding inflationary Universe the deviations in the spectrum were found BM1 ; MB2 , namely the (generalization of the) Corley/Jacobson dispersion relation (5) with $`b_m`$ negative. In this case, $`\omega _{_{\mathrm{𝐩𝐡𝐲𝐬}}}(k)`$ is given by

$$\omega _{_{\mathrm{𝐩𝐡𝐲𝐬}}}^2(k)=k^2|b_m|k^2\left(\frac{k^2}{k__\mathrm{C}^2}\right)^m,$$
(12)
where $`k__\mathrm{C}`$ is the cutoff physical wavenumber. The dispersion relation is represented in Fig. 1.

The behavior of the solutions of the equation of motion $`\mu ^{\prime \prime }+\mu (n_{\mathrm{eff}}^2a^{\prime \prime }/a)=0`$ is determined by the competition between the two terms $`n_{\mathrm{eff}}^2`$ and $`a^{\prime \prime }/a`$ given by

$$n_{\mathrm{eff}}^2(n,\eta )=n^2n^2\frac{|b_m|}{\mathrm{}__\mathrm{C}^{2m}}\left(\frac{n}{2\pi }\right)^{2m}\left[\mathrm{}_0\frac{\mathrm{}_0\mathrm{}_\mathrm{b}}{1+(\eta /\eta _0)^2}\right]^{2m},\frac{a^{\prime \prime }}{a}=\frac{2(\mathrm{}_0\mathrm{}_\mathrm{b})}{\eta _0^2\mathrm{}_\mathrm{b}}\frac{13(\eta /\eta _0)^2}{[1+(\eta /\eta _0)^2]^3},$$
(13)
where $`\mathrm{}__\mathrm{C}1/k__\mathrm{C}`$ is the cutoff length. These two terms are represented in Fig. 2. When $`n_{\mathrm{eff}}^2>a^{\prime \prime }/a`$, which is always the case when $`\eta /\eta _0\pm \mathrm{}`$, the WKB approximation is valid and the fundamental solutions can be written as

$$\mu \frac{1}{\sqrt{2n_{\mathrm{eff}}}}\mathrm{exp}\left[\pm i^\eta n_{\mathrm{eff}}(n,\tau )d\tau \right]_{\eta /\eta _0\pm \mathrm{}}\frac{1}{\sqrt{2n}}e^{\pm in\eta }.$$
(14)
In this case, the physical wavenumbers of the modes are such that they correspond to the linear part of the dispersion relation, see Fig. 1. In this situation, natural initial conditions can be chosen. Note that these initial conditions are, in a sense, even “more standard” than in the previous studies on trans-Planckian physics since usually the initial conditions are set in a region where the dispersion relation is not linear (but, as explained above, since the WKB approximation holds, meaningful initial conditions can nevertheless be considered). On the contrary, when $`n_{\mathrm{eff}}^2<a^{\prime \prime }/a`$ the WKB approximation is violated and two independent solutions are given by

$$\mu a(\eta ),\mu a(\eta )^\eta \frac{\mathrm{d}\tau }{a^2(\tau )}.$$
(15)
In this case, we are close to the cutoff scale and $`\omega _{_{\mathrm{𝐩𝐡𝐲𝐬}}}^2(k)`$ being close to zero, see Fig. 1, the $`a^{\prime \prime }/a`$ term dominates. This situation is very similar to the one discussed in Ref. Mersini  where this part of the dispersion relation has been named “the tail”. In fact, it is possible to quantify exactly the accuracy of the WKB approximation. Given an equation of the form $`\mu ^{\prime \prime }+\omega ^2\mu =0`$, the WKB approximation is valid if the quantity $`|Q/\omega ^2|1`$, where the quantity $`Q`$ is defined by the following expression

$$Q=\frac{3(\omega ^{})^2}{4\omega ^2}\frac{\omega ^{\prime \prime }}{2\omega }.$$
(16)
This standard criterion can be obtained in the following manner. The WKB solution, $`\mu _{\mathrm{wkb}}`$, satisfies the equation $`\mu _{\mathrm{wkb}}^{\prime \prime }+\mu _{\mathrm{wkb}}(\omega ^2Q)=0`$ exactly. Therefore, one has $`\mu \mu _{\mathrm{wkb}}`$ if $`|Q/\omega ^2|1`$. In the present context, $`\omega ^2`$ is of course equal to $`n_{\mathrm{eff}}^2a^{\prime \prime }/a`$. We have plotted the quantity $`|Q/\omega ^2|`$ in Fig. 3. When $`|\eta /\eta _0|1`$, $`|Q/\omega ^2|1`$ and the WKB approximation is satisfied, as expected from the previous discussion. The quantity $`|Q/\omega ^2|`$ increases when one approaches the bounce. At $`\eta =\eta _\mathrm{j}`$, where $`\omega (\eta _\mathrm{j})=0`$, there is a divergence. This point is a turning point. Then for $`|\eta /\eta _0|1`$, as can be seen in Fig. 3, $`|Q/\omega ^2|`$ remains large and the WKB approximation breaks down in agreement with the considerations above.

In a bouncing Universe, the modes start in the linear part of the dispersion relation, pass through the non-linear part and then come back in the standard region. However, things are not so simple because there exists a range of comoving wavenumbers such that the dispersion relation becomes complex. Although this case is a priori interesting, it clearly requires more speculative considerations, in particular the quantization in the presence of imaginary frequency modes. Since our goal in this paper is to exhibit a case where everything can be done in a standard manner, we will restrict ourselves to modes which never enter the region where the dispersion relation becomes complex. In the following, we determine the range of comoving wavenumbers we are interested in and for which we are going to calculate the power spectrum.

The maximum of the absolute value of $`a^{\prime \prime }/a`$ is located at $`\eta =0`$ and is given by $`(2/\eta _0^2)(\mathrm{}_0/\mathrm{}_\mathrm{b}1)`$. Therefore, if we restrict ourselves to modes such that

$$nn_\mathrm{b}\frac{\sqrt{2}}{\eta _0}\left(\frac{\mathrm{}_0}{\mathrm{}_\mathrm{b}}1\right)^{1/2},$$
(17)
then one is sure that, with an unmodified dispersion relation, the term $`a^{\prime \prime }/a`$ can always be neglected and that the initial spectrum is never changed. Clearly, this is not a physical restriction but it renders the comparison with the case with a modified dispersion relation easier. The minimum of the modified part of the dispersion relation is given by (from now on, we consider the case $`|b_m|=1`$ and $`m=1`$ since it is more convenient and does not restrict the physical content of the problem in any way)

$$n_{\mathrm{eff}}^2(n,\eta =0)=n^2\left[1\left(\frac{n}{n_{\mathrm{inf}}}\right)^2\right],n_{\mathrm{inf}}2\pi \frac{\mathrm{}_\mathrm{b}}{\mathrm{}__\mathrm{C}}.$$
(18)
Therefore, to maintain a real dispersion relation, we should only consider modes such that $`n<n_{\mathrm{inf}}`$. Of course, for consistency, the parameters of the model must be chosen such that $`n_\mathrm{b}<n_{\mathrm{inf}}`$. This is the case for the most natural choice, i.e., $`\mathrm{}_0\mathrm{}_\mathrm{b}\mathrm{}__\mathrm{C}`$. In the following, as already mentioned above, the time $`\eta _\mathrm{j}(n)`$ such that

$$n_{\mathrm{eff}}^2(n,\eta _\mathrm{j})=\frac{a^{\prime \prime }}{a}$$
(19)
will play a crucial role. For convenience we will only consider values of $`n`$ such that this time is determined only by the central peak of $`a^{\prime \prime }/a`$ and not by the two wings. In practice, this amounts to taking $`n_{\mathrm{eff}}^2(n,\mathrm{})>n_\mathrm{b}^2/4`$ since $`n_\mathrm{b}^2/4`$ is the maximum height of the wings. Then it is easy to show that $`n^2[n_{}^2,n_+^2]`$ where

$$n_\pm ^2=\frac{n_{\mathrm{sup}}^2}{2}\left(1\pm \sqrt{1\frac{n_\mathrm{b}^2}{n_{\mathrm{sup}}^2}}\right),$$
(20)
where $`n_{\mathrm{sup}}2\pi \mathrm{}_0/\mathrm{}__\mathrm{C}`$ is the largest value of $`n`$ such that $`n_{\mathrm{eff}}^2(n,\mathrm{})`$ remains real. Therefore, the range we are interested in is given by

$$\mathrm{max}(n_\mathrm{b}^2,n_{}^2)<n^2<\mathrm{min}(n_{\mathrm{inf}}^2,n_+^2).$$
(21)
In practice, we have $`n_{}^2(1/4)n_\mathrm{b}^2`$ and $`n_+^2n_{\mathrm{sup}}^2`$. This means that $`\mathrm{max}(n_\mathrm{b}^2,n_{}^2)=n_\mathrm{b}^2`$, $`\mathrm{min}(n_{\mathrm{inf}}^2,n_+^2)=n_{\mathrm{inf}}^2`$ and that the range reduces to

$$n_\mathrm{b}^2<n^2<n_{\mathrm{inf}}^2.$$
(22)
Having determined the relevant wavenumbers, we can now choose the initial condition and solve the equation of motion. In the first region where $`\eta <\eta _\mathrm{j}(n)`$, we only consider positive frequency modes and we have

$$\mu _\mathrm{I}(\eta )=\frac{1}{\sqrt{2n_{\mathrm{eff}}}}\mathrm{exp}\left[i_{\eta _\mathrm{i}}^\eta n_{\mathrm{eff}}(n,\tau )d\tau \right],$$
(23)
where $`\eta _\mathrm{i}`$ is an arbitrary initial time. In the second region, where $`\eta _\mathrm{j}(n)<\eta <\eta _\mathrm{j}(n)`$, the solution is given by

$$\mu _{\mathrm{II}}(\eta )=B_1a(\eta )+B_2a(\eta )_0^\eta \frac{\mathrm{d}\tau }{a^2(\tau )}.$$
(24)
The lower bound of the integral is a priori arbitrary. However, it is very convenient to take it equal to zero because in this case the second branch becomes odd whereas the first one (i.e. the scale factor) is even. Then, it is easy to show that

$$_0^\eta \frac{\mathrm{d}\tau }{a^2(\tau )}=\frac{\eta _0}{\mathrm{}_\mathrm{b}^2}\left\{\frac{1}{p^2}\frac{\eta }{\eta _0}+\frac{12p+p^2}{2p^2[1+p(\eta /\eta _0)^2]}\frac{\eta }{\eta _0}+\frac{p^2+2p3}{2p^{5/2}}\mathrm{arctan}\left(\sqrt{p}\frac{\eta }{\eta _0}\right)\right\},$$
(25)
where $`p\mathrm{}_0/\mathrm{}_\mathrm{b}`$. Finally, the solution in the third region where $`\eta >\eta _\mathrm{j}(n)`$ can be written as

$$\mu _{\mathrm{III}}(\eta )=\frac{C_1}{\sqrt{2n_{\mathrm{eff}}}}\mathrm{exp}\left[i_{\eta _\mathrm{i}}^\eta n_{\mathrm{eff}}(n,\tau )d\tau \right]+\frac{C_2}{\sqrt{2n_{\mathrm{eff}}}}\mathrm{exp}\left[+i_{\eta _\mathrm{i}}^\eta n_{\mathrm{eff}}(n,\tau )d\tau \right].$$
(26)
The goal is now to calculate the coefficients $`C_1`$ and $`C_2`$. Using the continuity of the mode function $`\mu `$ and of its derivative, we find

$`C_1(n)`$ $`=`$ $`{\displaystyle \frac{i}{2W[\eta _j(n)]}}{\displaystyle \frac{e^{i[\mathrm{\Omega }(\eta _\mathrm{j})\mathrm{\Omega }(\eta _\mathrm{j})]}}{\sqrt{n_{\mathrm{eff}}(n,\eta _\mathrm{j})n_{\mathrm{eff}}(n,\eta _\mathrm{j})}}}\left\{[g^{}+\alpha g](\eta _\mathrm{j})[f^{}+\overline{\alpha }f](\eta _\mathrm{j})+[f^{}+\alpha f](\eta _\mathrm{j})[g^{}+\overline{\alpha }g](\eta _\mathrm{j})\right\},`$ (27)
$`C_2(n)`$ $`=`$ $`{\displaystyle \frac{i}{2W[\eta _j(n)]}}{\displaystyle \frac{e^{i[\mathrm{\Omega }(\eta _\mathrm{j})+\mathrm{\Omega }(\eta _\mathrm{j})]}}{\sqrt{n_{\mathrm{eff}}(n,\eta _\mathrm{j})n_{\mathrm{eff}}(n,\eta _\mathrm{j})}}}\left\{[g^{}+\alpha g](\eta _\mathrm{j})[f^{}+\alpha f](\eta _\mathrm{j})+[f^{}+\alpha f](\eta _\mathrm{j})[g^{}+\alpha g](\eta _\mathrm{j})\right\},`$ (28)
where we have used the short-hand notation $`fa(\eta )`$ and $`ga(\eta )_0^\eta d\tau /a^2(\tau )`$. We have also utilized the following definitions: $`Wgf^{}g^{}f`$ is the Wronskian, $`\mathrm{\Omega }(\eta )_{\eta _\mathrm{i}}^\eta n_{\mathrm{eff}}(n,\tau )d\tau `$, and the quantity $`\alpha `$ is

$$\alpha \frac{n_{\mathrm{eff}}^{}}{2n_{\mathrm{eff}}}+in_{\mathrm{eff}}.$$
(29)
Using the fact that $`n_{\mathrm{eff}}`$ is even, it is easy to see that $`\alpha (n,\eta )=\overline{\alpha }(n,\eta )`$. Using the same property, one could also simplify the factor $`\sqrt{n_{\mathrm{eff}}(n,\eta _\mathrm{j})n_{\mathrm{eff}}(n,\eta _\mathrm{j})}=n_{\mathrm{eff}}(n,\eta _\mathrm{j})`$ in the above equation. The final result is given by Eqs. (27), (28) where all the functions are explicitly known except the function $`\eta _\mathrm{j}=\eta _\mathrm{j}(n)`$. The dependence on the dispersion relation of the final result is completely encoded in this function. The time $`\eta _\mathrm{j}(n)`$ is the solution of an algebraic equation that is not possible to solve explicitly. However, we can find an approximation if the function $`n_{\mathrm{eff}}(n,\eta )`$ is Taylor expanded around $`\eta =0`$. For the function $`a^{\prime \prime }/a`$ , we use a quadratic least square approximation which gives a better result. We obtain

$$n_{\mathrm{eff}}^2(n,\eta )n^2n^2\left(\frac{n}{n_{\mathrm{inf}}}\right)^2(1n_\mathrm{b}^2\eta ^2),\frac{a^{\prime \prime }}{a}n_\mathrm{b}^2\left[16\left(\frac{\eta }{\eta _0}\right)^2\right].$$
(30)
The two functions and their approximations are represented in Fig. 4. We see now that $`\eta _\mathrm{j}(n)`$ is the solution of a quadratic algebraic equation. Straightforward calculations give

$$\eta _\mathrm{j}^2(n)\frac{n^2}{n_\mathrm{b}^2}\left[\left(\frac{n}{n_{\mathrm{inf}}}\right)^2+\left(\frac{n_\mathrm{b}}{n}\right)^21\right]\left[n^2\left(\frac{n}{n_{\mathrm{inf}}}\right)^2+\frac{6}{\eta _0^2}\right]^1.$$
(31)
The function $`\eta _\mathrm{j}(n)`$ can also be found exactly by numerical calculations.

We have done this computation by means of a Fortran code. The comparison of the exact result with its approximation given in Eq. (31) is represented in Fig. 5. The function $`n_{\mathrm{eff}}(n,\eta =0)`$ is not monotonic and has a maximum around $`n45`$ for the values of the parameters considered in Fig. 5. For this value $`n_{\mathrm{eff}}(n,\eta )`$ is almost never smaller than $`a^{\prime \prime }/a`$. It is clear that around this value the approximation will be very good whereas for other values of $`n`$, especially for $`nn_\mathrm{b}`$, the approximation will be less good. This is due to the fact that when $`n`$ approaches $`n_\mathrm{b}`$, the asymptotic value of $`n_{\mathrm{eff}}(n,\eta )`$ decreases and the curve opens out at the intersection with $`a^{\prime \prime }/a`$. As a consequence, the quadratic approximation employed above breaks down. This is confirmed by the plots in Fig. 5. The approximation can be less good for other values of the parameters.

From the above analysis, it is clear that $`ϵ\eta _\mathrm{j}(n)/\eta _0`$ is a small number. We can therefore expand everything in terms of the parameter $`ϵ`$. This will allow us to get an analytical estimate of the spectrum. After tedious but straightforward computations, one obtains the following expressions for the coefficients $`C_1(n)`$ and $`C_2(n)`$

$`C_1(n)`$ $`=`$ $`1+ϵ{\displaystyle \frac{i\eta _0}{2n}}\left[2n^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{1/2}+n_\mathrm{b}^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{1/2}+n_\mathrm{b}^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{3/2}\right]+𝒪(ϵ^2),`$ (32)
$`C_2(n)`$ $`=`$ $`ϵ{\displaystyle \frac{i\eta _0}{2n}}e^{2i\mathrm{\Psi }_\mathrm{i}(n)}\left[2n^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{1/2}+n_\mathrm{b}^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{1/2}+n_\mathrm{b}^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{3/2}\right]+𝒪(ϵ^2),`$ (33)
where $`\mathrm{\Psi }_\mathrm{i}(n)_{\eta _\mathrm{i}}^0n_{\mathrm{eff}}(n,\tau )d\tau `$. Several comments are in order at this point. Firstly, the expressions of the two coefficients have exactly the expected form. When the parameter $`ϵ`$ goes to zero, $`C_1(n)`$ goes to one and $`C_2(n)`$ goes to zero, i.e., we recover an unmodified spectrum. This is in complete agreement with the fact that, when $`ϵ=\eta _\mathrm{j}(n)/\eta _00`$, the effective wavenumber $`n_{\mathrm{eff}}`$ never penetrates the region where the WKB approximation is violated. This is in this sense that the parameter $`ϵ`$ is the relevant parameter to expand the spectrum in since it is clearly directly linked to the violation of the WKB approximation. Therefore, in the limit $`ϵ0`$, one recovers the usual spectrum. Secondly, it is interesting to notice that the two terms proportional to $`ϵ`$ in $`C_1(n)`$ and $`C_2(n)`$ are the same up to a global phase. To compute the phase in Eqs. (27) and (28), we have written $`\mathrm{\Omega }(\eta )=\mathrm{\Psi }_\mathrm{i}(n)+_0^\eta n_{\mathrm{eff}}(n,\tau )d\tau `$. Since we only need to evaluate $`\mathrm{\Omega }(\pm \eta _\mathrm{j})`$, one can use the approximate equation for $`n_{\mathrm{eff}}`$, see Eq. (30), in the integral $`_0^{\eta _\mathrm{j}}n_{\mathrm{eff}}(n,\tau )d\tau `$. One finds

$$\mathrm{\Omega }(\pm \eta _\mathrm{j})=\mathrm{\Psi }_\mathrm{i}\pm ϵn\eta _0\left(1\frac{n^2}{n_{\mathrm{inf}}^2}\right)^{1/2}+𝒪(ϵ^3).$$
(34)
Let us remark that the term $`e^{i[\mathrm{\Omega }(\eta _\mathrm{j})+\mathrm{\Omega }(\eta _\mathrm{j})]}`$ in the expression of $`C_2(n)`$ does not contribute to order $`ϵ`$ whereas the factor $`e^{i[\mathrm{\Omega }(\eta _\mathrm{j})\mathrm{\Omega }(\eta _\mathrm{j})]}`$ gives a contribution at this order. This contribution is exactly such that the two terms proportional to $`ϵ`$ in Eqs. (32), (33) have the same absolute value. Thirdly, the factor $`1n^2/n_{\mathrm{inf}}^2`$ is always positive since we restrict ourselves to the regime where $`n^2<n_{\mathrm{inf}}^2`$. Therefore, the above expressions of $`C_1(n)`$ and $`C_2(n)`$ are always well-defined. Finally, interestingly enough, the two coefficients $`C_1(n)`$ and $`C_2(n)`$ have a similar structure to one of the corresponding coefficients found in Ref. LLMU  where another example of a modified spectrum has been exhibited.

We are now in a position where we can calculate the spectrum explicitly. It is defined by the following equation

$$n^3P(n)n^3\left|\frac{\mu }{a}\right|^2.$$
(35)
We will evaluate the spectrum in the region where $`\eta /\eta _01`$. In this case, $`a(\eta )\mathrm{}_0`$ and $`n_{\mathrm{eff}}n`$. Expanding everything in terms of the small parameter $`ϵ`$, one obtains

$`n^3P(n)`$ $``$ $`{\displaystyle \frac{n^2}{2\mathrm{}_0^2}}+ϵ{\displaystyle \frac{\eta _0}{2\mathrm{}_0^2}}n\left[2n^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{1/2}+n_\mathrm{b}^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{1/2}+n_\mathrm{b}^2\left(1{\displaystyle \frac{n^2}{n_{\mathrm{inf}}^2}}\right)^{3/2}\right]`$ (36)
$`\times \mathrm{sin}\left[2\mathrm{\Phi }(n)2\mathrm{\Psi }_\mathrm{i}(n)+2n(\eta \eta _+)\right]+𝒪(ϵ^2),`$
where $`\eta _+>\eta _\mathrm{j}`$ is the time at which $`n_{\mathrm{eff}}n`$ and the phase $`\mathrm{\Phi }(n)`$ is defined by $`\mathrm{\Phi }(n)_{\eta _\mathrm{i}}^{\eta _+}d\tau n_{\mathrm{eff}}`$. Again the spectrum has the expected form. The leading order is the vacuum spectrum with spectral index $`3`$ and the next-to-leading order (proportional to $`ϵ`$) has the form of a complicated function of $`n`$ times superimposed oscillations. We have plotted this complicated function (without the superimposed oscillations) in Fig. 6.

We can check that the correction remains small and that, therefore, the approximation employed in this article is consistent. On the other hand, when the wavenumber $`n`$ approaches $`n_{\mathrm{inf}}`$, the approximation scheme breaks down but, precisely in this limit, the whole problem becomes pointless in agreement with the considerations above. We also see that the deviation is minimum around $`n45`$ (for these values of the parameters) because, as already mentioned above, in this case the term $`a^{\prime \prime }/a`$ almost never dominates. Although the present example is only a toy model, it is quite interesting to see that it bears some close resemblance with the other example of a modified spectrum found in Ref. LLMU . This suggests that the features previously exhibited may be valid in general.

We have thus reached our main goal: find an example where the initial conditions can be fixed naturally, where the frequency never becomes complex and where the final spectrum is modified. Note that in this example, the change in the spectrum comes about from an interplay between the modified dispersion relation factor $`n_{\mathrm{eff}}^2`$ and the factor $`a^{\prime \prime }/a`$ which is responsible for Parker particle production Parker . For an unmodified linear dispersion relation, the $`a^{\prime \prime }/a`$ term is always negligible and there is no particle production. However, for our modified dispersion relation, for a range of modes the Parker particle production term becomes important and leads to non-adiabatic evolution. The length of time during which the $`a^{\prime \prime }/a`$ term dominates depends on the specifics of the dispersion relation, and hence the final spectrum will depend on these specifics.

## IV Discussion

In this paper we have studied the dependence of the spectrum of a free scalar field in a bouncing Universe on trans-Planckian effects introduced via a modification of the free field dispersion relation. We have found that both the amplitude and the slope of the spectrum depend on the dispersion relation, assuming that the dispersion relation leads to non-adiabatic evolution of the mode functions on trans-Planckian physical length scales. Such non-adiabatic evolution is possible without requiring the effective frequency to become imaginary.

This result supports our earlier work which shows that the spectrum of cosmological fluctuations in inflationary cosmology can depend sensitively on trans-Planckian physics BM1 ; MB2 . Note that results supporting BM1 ; MB2  were recently also obtained by EGKS  (see also Ref. KN ) in the context of a mode equation modified from the usual ones by taking into account Kempf  effects coming from a short distance cutoff. For a recent paper on the effects of short-distance physics on the consistency relation for scalar and tensor fluctuations in inflationary cosmology see Ref. Hui . The calculation we presented here is free of two of the possible objections against the earlier work. The first objection was that the use of the minimum energy density state as initial state is not well defined on trans-Planckian scales if the dispersion relation differs dramatically from the linear one (which it has to in order to get non-adiabatic evolution). In our present work, the initial conditions are set in the low curvature region and on length scales larger than the cutoff length but smaller than the Hubble radius, where the choice of initial vacuum state is well defined. The second objection concerned the use of dispersion relations for which the effective frequency is imaginary in some time interval. No imaginary frequencies are used in the present model.

The results obtained in this paper apply also both to bouncing Universe backgrounds described by the two other models mentioned in the Introduction. In the case of the first model (where the scale factor is given by an hyperbolic cosine), similar results would hold if we would match the local Minkowski vacuum (WKB vacuum) at some initial time and express the results in terms of the WKB vacuum state at the corresponding post-bounce time. However, in that case, the initial conditions are not easy to justify since they have to be set on length scales larger than the Hubble radius. Obviously, we could consider that model and let the scale factor make a further transition to an asymptotically radiation- dominated or asymptotically flat stage at very large initial and final times. In this case, it is very likely that results similar to the ones we obtained here in model (1) would be obtained.

Acknowledgments

We would like to thank M. Lemoine for comments and careful reading of the manuscript. We acknowledge support from the BROWN-CNRS University Accord which made possible the visit of J. M. to Brown during which most of the work on this project was done, and we are grateful to Herb Fried for his efforts to secure this Accord. S. E. J. acknowledges financial support from CNPq. R. B. wishes to thank Bill Unruh for hospitality at the University of British Columbia during the time when this work was completed, and for many stimulating discussions. J. M. thanks the High Energy Group of Brown University for warm hospitality. The research was supported in part by the U.S. Department of Energy under Contract DE-FG02-91ER40688, TASK A.