{'fraction_non_alphanumeric': 0.030435395237424265, 'fraction_numerical': 0.008383824151049738, 'mean_word_length': 4.355005030181086, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Instability of the Lee-Wick Bounce

It was recently realized that a model constructed from a Lee-Wick type scalar field theory yields, at the level of homogeneous and isotropic background cosmology, a bouncing cosmology. However, bouncing cosmologies induced by pressure-less matter are in general unstable to the addition of relativistic matter (i.e. radiation). Here we study the possibility of obtaining a bouncing cosmology if we add radiation coupled to the Lee-Wick scalar field. This coupling in principle would allow the energy to flow from radiation to matter, thus providing a drain for the radiation energy. However, we find that it takes an extremely unlikely fine tuning of the initial phases of the field configurations for a sufficient amount of radiative energy to flow into matter. For general initial conditions, the evolution leads to a singularity rather than a smooth bounce.


I. INTRODUCTION
Both Standard [2] and Inflationary Cosmology [3] suffer from the initial singularity problem and hence cannot yield complete descriptions of the very early universe. If one were able to construct a non-singular bouncing cosmology, this problem would obviously disappear. However, in order to have a chance to obtain such a nonsingular cosmology, one must either go beyond Einstein gravity as a theory of space-time (see e.g. [4] for an early construction), or else one must make use of matter which violates the "null energy condition" (see [5] for a review of both types of approaches).
Interest in non-singular bouncing cosmologies has increased with the realization that they can lead to alternatives to inflationary cosmology as a theory for the origin of structure in the universe. A specific scenario which can arise at the level of homogeneous and isotropic cosmology is the "matter bounce" paradigm which is based on the realization [6,7] that vacuum fluctuations which exit the Hubble radius during a matter-dominated contracting phase evolve into a scale-invariant spectrum of curvature perturbations on super-Hubble scales before the bounce. The key point is that the curvature fluctuation variable ζ grows on super-Hubble scales in a contracting phase, whereas it is constant on these large scales in an expanding phase. Since long wavelength modes exit the Hubble radius earlier than short wavelength ones, they grow for a longer period of time. This provides a mechanism for reddening the initial vacuum spectrum. It turns out that a matter dominated contracting phase provides the specific boost in the power of long wavelength modes which is required in order to transform a vacuum spectrum into a scale-invariant one. Studies in the case of various non-singular bounce models [8] have shown that on wavelengths long compared to the duration of the bounce phase, the spectrum of fluctuations is virtually unchanged during the bounce. Thus, a scaleinvariant spectrum of curvature fluctuations survives on super-Hubble scales at late times.
Provided that the bounce can occur at energy scales much below the Planck scale, non-singular cosmologies solve a key conceptual problem from which inflationary cosmology suffers, namely the "Trans-Planckian" problem for fluctuations [9,10]: If the period of inflationary expansion of space lasts for more than 70H −1 , where H is the Hubble expansion rate during inflation (in order to solve the key cosmological mysteries it was designed to explain, inflation has to last at least 50H −1 ), then the physical wavelengths of even the largest-scale fluctuation modes we see today will be even smaller than the Planck length at the beginning of inflation and thus in the "zone of ignorance" where the physics on which inflation and the theory of cosmological perturbations are based, namely Einstein gravity coupled to semiclassical field theory matter, will break down. In contrast, in a non-singular bouncing cosmology the wavelength of modes which are currently probed by cosmological observations is never much smaller than 1mm (the physical wavelength of the mode which corresponds to our current Hubble radius evaluated when the temperature of the universe was 10 16 GeV) and hence many orders of magnitude larger than the Planck length. Thus, the fluctuations never enter the "trans-Planckian zone of ignorance" of sub-Planck-length wavelengths.
Possibly the simplest realization of the matter bounce scenario is the "quintom bounce" model [11] and is obtained by considering the matter sector to contain two scalar fields, one of them (the "ghost field") having the "wrong" sign of the kinetic action. The potential of the ghost scalar field also has the opposite sign to that of regular scalar fields such that in the absence of interactions, the ghost field has a classically stable minimum. As has been noticed in [1], such a quintom bounce model also arises from the scalar field sector of the "Lee-Wick" (LW) Lagrangian [12] which contains higher derivatives terms.
The quintom and Lee-Wick bouncing cosmologies are obtained in the following way [1,11]: We begin in the contracting phase with both the regular and the ghost scalar field oscillating homogeneously in space about their respective vacua. We assume that the energy density is dominated by the regular matter field, and that hence the total energy density is positive. Once the amplitude of the regular scalar field exceeds the Planck scale, the field oscillations will freeze out and a slowclimb phase will begin during which the energy density of the field only grows slowly (this is the time reverse of the slow-roll phase in scalar field-driven inflation). However, the ghost field continues to oscillate and its energy density (which is negative) continues to grow in absolute value. Hence, the total energy density drops to zero, at which point the bounce occurs, as has been studied both analytically and numerically in the abovementioned works. Note that the energy density in this bounce model scales as matter until the regular scalar field freezes out.
A major problem of bouncing cosmologies realized with matter which scales as a −3 as a function of the scale factor a(t) is the potential instability of the homogeneous and isotropic background against the effects of radiation (which scales as a −4 and anisotropic stress which scales as a −6 [16]. If we simply add a non-interacting radiation component to the two scalar field system, then unless the initial energy density in radiation is tuned to be extremely small, then the radiation component will become dominant long before the bounce can arise, and will prevent the energy density in the ghost field from ever being able to become important, resulting in a Big Crunch singularity. Similarly, unless the initial energy density in anisotropic stress is very small, it will come to dominate the energy density of the universe long before the bounce is expected. The anisotropies will destabilize the homogeneous background cosmology, and will prevent a bounce. Note that at the quantum level, there is an additional severe problem for bounce models obtained with matter fields with ghost-like kinetic terms, namely the quantum instability of the vacuum (see e.g. [14]).
In this paper we will focus on the radiation instability problem. For the purpose of this discussion we will simply assume that anisotropic stress is absent. In a recent paper, two of us studied the possibility that a bounce could arise if radiation is supplemented with Lee-Wick radiation [15]. However, we showed that this hope is not realized: the addition of Lee-Wick radiation does not prevent the Big Crunch singularity from occurring. In the presence of radiation, the only hope to obtain a bounce is to introduce a coupling between radiation and ghost scalar field matter which could effectively drain energy density from the radiation field and prevent the energy density of radiation from becoming dominant. Here we study this possibility. However, at least for the specific Lagrangian which we consider, we find that a bounce only emerges for highly fine-tuned phases of the fields and their velocities in the initial conditions. The paper is organized as follows: In Section II, we introduce the model we study, namely the scalar field sector of Lee-Wick theory coupled to radiation, and write down the general equations of motion. In Section III we set up the equations of motion linearized about the bounce background, treating the entire radiation field as an inhomogeneous fluctuation. In particular, we study the different terms which contribute to the energy-momentum tensor and identify those which could assist in obtaining a non-singular bounce. In Section IV we study the solutions of the perturbed equations of motion, and in Section V we analyze the evolution of the different terms in the energy-momentum tensor, identifying the conditions which would be required in order to obtain a non-singular bounce. We have also evolved the general equations of motion for the two inhomogeneous scalar field configurations and the classical inhomogeneous radiation field in the homogeneous background cosmology. Section VI summarizes some of the numerical results. Both the analytical and numerical results confirm that we need unnatural fine-tuning of the initial conditions in order to obtain a non-singular bounce. In the final section we offer some conclusions and discussion.

II. THE MODEL
The Lee-Wick scalar field model coupled to electromagnetic radiation is given by the following Lagrangian: where m is the mass of the scalar field φ, and V (φ) is its potential. Here we adopt the convention that where a(t) is the scale factor of the universe. Since it is a higher derivative Lagrangian in φ, the scalar field sector contains an extra degree of freedom with the "wrong" sign kinetic term and with a mass set by the scale M . We choose m ≤ M ≤ m P l , where m P l is the Planck mass, since we want the regular scalar field to dominate at low energies, but at the same time we do not want to worry about quantum gravity effects. The second line of the Lagrangian (1) contains the kinetic term of the radiation as well as the coupling term, where we assumed both for the sake of generality and because of foresight that the radiation field couples not only to the scalar field φ itself, but also to the higher derivative term. The electromagnetic tensor, F µν , is related to the radiation field A µ through the usual definition where ∇ µ is the covariant derivative.
It is convenient to extract the extra degree of freedom as a separate scalar field. To do this, we use the field redefinitions The Lagrangian (1) then takes on a simpler form: where we have chosen the potential to be zero. In this new form, the Lagrangian describes two massive scalar fields with one of them (i.e., φ 2 ) behaving like a"ghost", and both of them coupled to the radiation field.
The coupling term f (φ 1 , φ 2 , F µν F µν ) should in principle be arbitrary, however, in this paper we will take a specific form for convenience. The form will be: where c and d are coupling constants which have mass dimension −2. The interaction terms are nonrenormalizable. To make sure that such terms could be thought of as arising from an effective field theory which is consistent at the bounce, we must make sure that the coefficients are chosen such that the contribution of the interaction term to the Lagrangian density is smaller than that of the other terms. This must be true even at energy densities at which the bounce occurs in the pure scalar field model. It is easy to see that this condition will be satisfied if the coefficients c and d are both of the order m −2 pl . It is the purpose of this paper to study the effects which these coupling terms have on the dynamics of the system. We know that in the absence of coupling, i.e. when c = d = 0, a bounce will only occur if the initial radiation energy density is tuned to a very small value compared to the scalar field energy density. This is because the positive definite energy density of radiation will scale as a −4 which is faster than that of the scalar fields, in particular the ghost scalar field. Generically, it will dominate the energy of the universe after some amount of contraction, it will prevent the ghost scalar field energy density from catching up and will thus prevent a bounce, leading to a Big Crunch singularity instead. With non-vanishing values of c and d, however, the scalars are in principle able to drain energy from the radiation.
From the Lagrangian (5), one can obtain the stressenergy tensor T µν by varying the action with respect to the metric g µν . In the hydrodynamical limit, we can take T µν to be of the form of diag{ρ, a 2 (t)p 1 , a 2 (t)p 2 , a 2 (t)p 3 } where ρ and p are energy density and pressure, respectively. For now, we consider the general form which is: Since we will be studying the contribution of plane wave perturbations of the scalar fields and we will treat radiation as a superposition of waves, we kept the spacederivative terms.
By varying the Lagrangian with respect to the matter fields φ 1 , φ 2 and A µ , we also get the equations of motion for all three fields: which will be analyzed in detail in the rest of the paper.

III. DYNAMICS
Since the equations of motion are nonlinear, we cannot work in Fourier space, and use plane wave solutions. However, we are interested in how initially small amounts of radiation build up and possibly transfer their energy to scalar field fluctuations. We treat radiation as a superposition of fluctuations. Therefore it makes sense to linearize our equations about the homogeneous scalar field background. Thus, we make the following ansatz for the scalar fields: 2 (t) , (12) where the expansion parameter ǫ is taken to be much smaller than 1 [17]. The first term on the right hand side of each line, i.e. φ To simplify the analysis, we describe radiation in terms of plane waves in a fixed direction (which we take to be the z direction). Without loss of generality we can restrict attention to one polarization mode which we take to be the electric field in the x direction and the magnetic field in the y direction. In this case, the only nonzero components of the field strength tensor are F 01 and F 13 . Using the temporal gauge where A 0 = 0, we find that only the first component of the gauge field is nonzero. For a single wavelength fluctuation we can make the ansatz or, equivalently, Since in the linearized equations of motion the Fourier modes are independent, we can consider φ (1) 1 and φ (1) 2 also to be plane waves propagating in z direction, so they depend only on z and t. With Eqs. (11)(12)(13)(14) in hand, we can write down the energy densities of the various fields at each order in perturbation theory.

A. The stress-energy tensor
First of all, we insert the above perturbative ansatz for the fields into the stress-energy tensor of the system. From the general expression (7) for T µν we get: The 00 component of Eq. (15) denotes the energy density of the system: so at each level in perturbation theory we have: 2 ) 2 ) We can similarly obtain the pressure of the system from the ii components of Eq. (15). Note that due to the anisotropy in T µν caused by the gauge field as well as by the anisotropic fluctuations of the scalar fields, the pressures in the three directions are no longer identical. The pressure in each direction can be written as: with no summation over the index i. From this formula, we can see that at both zero-th and first order, the pressure is isotropic: while the second order pressure for each direction gives where i = 1, 2, 3. We can thus obtain every component of p (2) i : 2 ) 2 ) and the average is: 2 ) From the above, we can also see that in order to analyze the behavior of the energy density up to second order, we need to know the evolution of scalar fields up to second order as well as that of the gauge field up to first order, while the behavior of the gauge field to second order is not required.
It will be useful in the following to separate the contributions to the energy density and pressure in a different way, namely, i) the contribution from the background homogeneous part of the scalar fields, ii) that of the scalar field perturbations (in slight abuse of notation we call this the "inhomogeneous" term), iii) the contribution of the gauge field, and iv) the contribution of the coupling term, where in the last equation we define Φ to be the quadratic combination of the two fields: From the above we can deduce the equation of state parameter for each part: From the equations above, we see that for positive values of the constants c and d, the coupling of the scalar field with the gauge field will give rise to a contribution ρ c to the energy density which has the same equation of state but opposite sign to that of the gauge field. Therefore the coupling can help drain energy from the gauge field. It is because of this mechanism that we might hope to achieve a cosmological bounce in the presence of radiation. A first indication on whether a bounce might occur can be obtained by considering the scaling of each contribution to the energy density as a function of the scale factor a(t). To find these scalings, we need the time dependence of the linear and quadratic contributions to each field. Therefore, we need to solve the matter field equations of motion. In the following subsection, we present the equations for the fields at each order, while the solutions and detailed analysis will be performed in the next sections.

B. Equations of Motion
Keeping in mind the ansätze for A µ , φ 1 and φ 2 , their equations of motion at each order can be obtained from (8), (9) and (10): a) At zero-th order: Note that there is no equation at this order for A µ because it is of first order in ǫ. b) At first order: Making use of Eqs. (13) and (14), the equation for the gauge field can also be rewritten as: c) At second order: Here, pointed parentheses indicate spatial averaging (since we are only focusing on the zero mode of the second order field fluctuations). We also neglected the effect of Hubble friction since it does not give an important contribution for the second order fluctuations.

IV. THE GENERAL SOLUTION
In this section we will solve the equations of motion (39), (40) and (42) to see if and how a bounce will happen.
It is usually useful to perform the analysis in the conformal frame where the conformal time η ≡ a −1 (t)dt is used rather than the cosmic time. Additionally, to extract the dependence on the scale factor, it is convenient to use the following two variables: Hereafter, we will use u (i) j (i = 0, 1, 2, j = 1, 2) to denote the i-th order perturbation of the j-th scalar field. Moreover, for simplicity but without loosing generality, we can parameterize the scale factor a(t) as where a 0 and w are the initial value of the scale factor and the equation of state of the universe, respectively. This is a self-consistent assumption when w is nearly a constant. The evolution of w in our case will be shown numerically in Section VI.
A. Solutions for φ Using the parametrization (44), the equations of motion at zeroth order of the two scalar fields become: where a prime denotes the derivative with respect to conformal time η. Their solutions are: where H ± 1−3p 2 represents the (± 1−3p 2 )-th order Hankel function. Far away or close to the bounce, i.e. for a|η| ≫ m −1 , M −1 and a|η| ≪ m −1 , M −1 , respectively, the approximate solutions are: 1)Oscillations for large values of the scale factor a|η| ≫ m −1 , M −1 : (where θ 1 and θ 2 are phases set by the initial conditions).
In terms of the non-rescaled fields φ i one obtains damped (or anti-damped) oscillations (depending on whether we are in an expanding or a contracting period) 2)"Frozen" evolution for small values of the scale factor a|η| ≪ m −1 , M −1 : from which it follows that the non-rescaled fields φ i evolve as from which we can see that the last term of φ (0) i is a constant mode while the first term is a varying one. Depending on the value of p (or equivalently w) the varying mode could be growing (for p > 1/3 or −1 < w < 1), in which case it becomes dominant, or decaying (for p < 1/3 or w > 1/w < −1), in which case it becomes subdominant. We can usually neglect the decaying part of the fields. Following the steps performed in the last subsection, we can also get the solutions for the first order components of the scalar fields. Using the equations (40) for the first order perturbations we obtain the following equations of motion for u Depending on the value of k, we obtain different approximation solutions. For wavenumbers large compared both to the Hubble radius and to the mass term, we obtain oscillatory solutions with fixed amplitude.
Considering now modes which are still sub-Hubble (i.e. k|η| > 1) but for which the mass term dominates over the contribution of the field tension (i.e. the term involving k), we can neglect both the k 2 term and the term involving p(2p−1) (1−p) 2 η 2 . The simplified equation for these modes is: whose solutions are: For modes outside the Hubble radius (kη ≪ 1), we have: which have the same form as Eq. (46) so their solution will be the same as given in Eqs. (53) and (54).
We have thus seen that the first order solutions for the scalar fields scale the same way with |η| as the zeroth order solution. This is because in the small |η| region where the a ′′ (t)/a(t) term dominates over the other ones, the equations for first order and zero-th order modes are almost the same. Thus, unless the energy density in the u (1) i modes dominates at the initial time, it will never dominate over the background contribution from the u (0) i terms. Thus we can conclude that the first order fluctuations of scalar fields will not prevent the bounce.

C. Solution for the gauge field γ
In this section, we will analyze the gauge field γ which is also considered to be of first order. The equation (41) can directly be transformed to conformal frame as: Since the coefficients c and d are small, we can take the last term to be a source term. In a first order Born approximation, we can write the total solution as where γ 0 is the solution for the homogeneous equation obtained by setting c = d = 0, while δγ is the leading correction term obtained by inserting γ 0 into the source term (the last term in (62)). The zero-th order (homogeneous) equation is easily solved and gives For the first order equation, it is convenient to define so that the equation becomes where we neglected the small term P (η)δγ ′ . Inserting the solution of γ 0 (64), we get the following equation for δγ: We are interested in the scaling of δγ as a function of time. For this purpose, we need to work out the scaling in time of the source term in (67). Since the solutions for φ For times obeying a|η| ≫ m −1 , M −1 , then by differentiating (51) and (52) with respect to η we have: Note that |η| 1 1−p ∼ t is a decaying mode in the contracting phase and thus the last terms inside the square brackets in the above formulae can be neglected compared to the first ones. Since then combining all these results we get: For a|η| ≪ m −1 , M −1 , then differentiating (55) and (56) with respect to η we obtain when p > 1/3 and when p < 1/3. Then we can solve Equation (67) to get: In the above expressions for δγ, C 1 , C 2 and C 3 are complicated prefactors in front of the η-dependent terms.
In summary, we see that the interactions give only a subleading correction δγ to γ 0 . 2 . If we only consider the energy density up to second order, these second order field perturbations give a contribution through their coupling to the background fields. In the following we find the solutions of (42) and study the effects of the induced terms in the stress-energy tensor on a possible bounce.

D. Solutions for φ
Given the solution for the gauge field γ obtained in the last subsection, it is easy to rewrite Eqs. (42) as: where we made use of the fact that Equation (76) has the same form as the zero-th order equation but with a small source term generated by the interaction with the gauge field. This equation can be solved using the Born approximation (details are given in the Appendix). The general solution is the sum of the general solution of the homogeneous solution plus the solution including the source which has vanishing initial data. The inhomogeneous term is suppressed by the coupling constants c and d compared to the homogeneous solution, but, as shown in the Appendix, it scales as a high power of η −1 . Via the coupling to the background scalar fields, the above second order terms enter into the expression for the energy density to second order. The signs of the corresponding terms in the energy density are indefinite in the sense that they depend on the phases of the initial field configurations. Since it is these terms that dominate the energy density near the bounce, we find that whether a bounce occurs or not depends sensitively on the phases in the initial conditions, and that in fact in the case of many plane wave modes initially excited, a bounce requires very special phase correlations.

V. EVOLUTION OF THE COMPONENTS OF THE ENERGY DENSITY
In the previous section we have solved all of the field equations up to second order in the amplitude of the fluc-tuations. We have found the scaling in time of each field at each order. Now we are ready to look at how all of the terms in the expression for the energy density ρ (0) , ρ (1) and ρ (2) at various orders in perturbation theory (namely, Eqs. (17)-(19)) scale in time. This analysis is straightforward but very important if we are to determine whether a bounce is possible, since in four space-time dimensional classical Einstein Gravity with flat spatial sections a bounce can only happen when the negative terms in the energy density catch up to the positive contributions [11].
In the following we give a table of how each term contained in ρ scales with time as the background cosmology bounce point (the bounce which is achieved in the absence of radiation and scalar field inhomogeneities) is approached. We will identify the terms which dominate in this limit. This will give us a good indication under which conditions a bounce can occur. The tables are structured as follows: the first line "Terms", indicates which term we are considering, the next set of lines "Behavior" gives the scaling in time in the various limits and in the two relevant ranges of the parameter p which indicates the equation of state, and the last line gives the sign with which the term contributes to the energy density. Note that we focus on the growing mode solution to each field (which is constant for small η in the case p < 1/3). We give separate tables for terms of zero-th, first and second order in ǫ. a) For terms contained in ρ (0) :

Sign Positive Definite
Indefinite (Depending only on c and d) Indefinite Indefinite These terms to not vanish upon spatial averaging. Note that we have expressed the time dependence in terms of the dependence on the scale factor a(t). At this stage, we only need to focus on the exponent of the power-law scaling. The more negative the power is, the more rapidly the term grows in a contracting phase (since a(t) is decreasing with time).
As mentioned earlier, the conditions for a bounce to occur in four space-time dimensional classical Einstein gravity with flat spatial sections is that the total energy density reaches zero during the contracting phase. Thus, there needs to be a negative definite term which starts out small but grows faster than the positive definite terms due to the regular scalar field and regular radiation. In the absence of radiation and scalar field inhomogeneities, it is the contribution to the energy density of the ghost field φ 2 which plays this role.
From the table we see that there are three kinds of terms: positive definite, negative definite and indefinite ones. The first set contains the kinetic and potential terms of the normal scalar as well as the free energy density of the gauge field, the second set is made up of the kinetic and potential terms of the ghost scalar, while the third set contains terms which arise due to the coupling terms between scalars or between scalars and gauge fields.
Looking first at the terms which are independent of the coupling term between the fields, we see from the first line of the "Behavior" set of lines in the third table that, indeed, in the presence of radiation the energy density in radiation grows faster than that in the two scalar fields, thus preventing a bounce. In the presence of coupling between the fields, however, there are terms which scale with a larger negative power of a(t). The signs of some of them, however, depend on the initial phases for the linear fields γ, φ Note that the signs of the scalar coupling terms are determined by the evolution of each field and thus are hard to be identified in a general analysis. The same is true for the gauge coupling terms (the last two in the third table). However, the coupling terms between the scalar fields and the gauge field (the third to last in the third table) can be made negative/positive definite easily by setting the signs of the coefficients c and d to be both positive/negative.
It is reasonable to assume that the contracting phase begins with the regular scalar field dominating the energy density, and that the contribution of the Lee-Wick scalar is much smaller. For single Fourier mode initial conditions of the radiation field, this can be achieved with the appropriate choice of the initial phase (see Example 1 in the following section containing our numerical results). However, for multiple initial radiation Fourier modes excited any initial phase difference between the modes will produce a contribution with the wrong sign and will thus prevent a bounce (see Example 3 in the following section). In the presence of an infinite set of modes, the phase correlations required to obtain a bounce thus appear to have negligible measure in initial condition space. Thus, even in the presence of coupling between scalar fields and radiation, the Lee-Wick bounce is unstable. The bounce, if it exists, will happen at a time which can be chosen to be t = 0. Its duration (the time interval lasting from the time the Hubble radius stops decreasing in the contracting phase until when it starts expanding in the post-bounce phase) will be denoted by ∆t. Since the various components of the energy density scale with different powers of a(t), it is clear that the duration of the bounce will be shorter or equal to the Hubble radius H −1 max (which gives the time scale on which the ratios of energy densities in different components change) at the beginning of the bounce phase. For the background bounce model, we have H max ∼ m.
There are two kinds of bounces according to the duration of the bounce phase i) If the period ∆t ≃ m −1 , the bounce will go from the time to the time with a low speed. We call this a "slow bounce". In this case, the universe will enter the bounce period at the critical time t c ∼ m, and only the a|η| ≫ m −1 approximate solutions of the previous tables will be applicable and not ones for the interval a|η| ≪ m −1 . ii) If the period ∆t ≪ m −1 , the bounce will happen in a very short time with very fast speed. This can be called the "fast bounce". In this case, the universe evolves from the far past (−t i with |t i | ≫ 1) to t = 0, passing through the point t c ∼ m, then entering into the region |t| m before finally reaching the bounce point. In this case, both of the two approximate solutions of the field evolution will be applied. Let us now consider the necessary conditions for a bounce (as we have indicated above and will see from the numerical analysis, these conditions are not sufficient -in addition to the conditions which follow, appropriate correlations in the initial phases are required). We start in the region of time a|η| ≫ m −1 . We study the conditions required to have the terms that might give a bounce grow relative to the other terms during this phase. If the conditions are not satisfied, or the bounce does not happen even if the conditions are satisfied, then a bounce may still occure in the a|η| ≪ m −1 region. The conditions for the terms in the energy which could compensate the positive radiation contribution to become dominant are then studied. If these conditions are not satisfied, either, then a bounce is impossible.
A necessary condition for a bounce to be possible requires the growth rate of one of the indefinite sign terms in the third table above exceed all that of all of the positive definite terms. In the a|η| ≫ m −1 region, this requires 8 − 1/p > 5, which equivalently constrains the equation of state parameter w to be in the range w < 1. If this condition is satisfied in this region, then a slow bounce may happen depending on the choice of the initial phases.
If the condition is not satisfied in the a|η| ≫ m −1 region, the universe may evolve into the a|η| ≪ m −1 region, in which the evolution of the fields are different, and new constraints on p and w will arise if a bounce is to be possible. Following the above logic, we find that the conditions under which a bounce might happen are much looser, namely w > −7/6.
To summarize this section: we have identified necessary conditions for a bounce to occur. Whether one actually does occur even if the conditions are satisfied depends on the initial phases of the fields. This must be studied numerically. In the following section we will give one example of specially chosen phases for which a bounce is possible. However, when we look at a more general choice of phases, the bounce will not occur.

VI. NUMERICAL RESULTS
In order to support the analysis in the last section, we performed numerical calculations. Such numerical work is necessary because our analytical analysis is only approximate. In particular, we worked in perturbation theory up to order second order in ǫ. In addition, even in cases where our analytical analysis would indicate the possibility of a bounce, the perturbative analysis will break down near the bounce point, and there is no assurance that the trends seen in the perturbative analysis will persist.
We have numerically solved the full nonlinear equa-tions of motion for the matter fields in the presence of a homogeneous expanding background cosmology. The homogeneous cosmology is obtained numerically by solving the first Friedmann equation where G is Newton's gravitational constant (related to the Planck mass used earlier), and ρ is the total energy density, averaged over space. Figures 1-6 are two groups of numerical results with different parameters. In both cases we choose the initial energy density of the gauge fields to be larger than that of the Lee-Wick scalar, but less than that of the normal scalar. These initial conditions correspond to the situation we are interested in, namely starting in a matterdominated contracting phase in the presence of some radiation which is sub-dominant. Fig. 1, 2 and 3 show an example with parameters c > 0 and d > 0. We choose initial conditions in which a single Fourier mode fluctuation is excited, and in which the phases are chosen as indicated in the figure caption. For these initial phases, we obtain a bounce. In Fig. 1, we see that the equation of state w begins with a value slightly larger than 0, and then evolves to some nearly fixed value. For the case of our initial condition choice, it appears to be w ≃ −0.6, in the region where the bounce is allowed to happen. At the bounce point, the equation of state will drop to −∞, while after the bounce, the equation of state will rise again to w ≃ 0.6. Fig. 2 is the plot of the scale factor in this case which shows explicitly the occurrence of the bounce. Fig. 3 gives a comparison of the energy densities of some components during the process. Initially, we set the energy density of the gauge field γ to be between the normal scalar and Lee-Wick scalar. When the evolution of the universe enters into a region with nearly constant w, the gauge-coupling component of energy density ρ c will grow very fast. It is negative and thus enables the negative part of the energy density to catch up with the positive one, thus allowing the bounce to happen. For the inhomogeneous fluctuation, we choose the wavenumber to be k ≃ 0.01hM pc −1 which corresponds to a scale which is observable by CMB and LSS experiments.
Figs. 4, 5 and 6 give the corresponding results in the case when we choose c > 0 while d < 0 (with all initial conditions identical). This case seems dangerous because the contribution of the Lee-Wick scalar to the fluctuation terms could lead to an instability. However, as we have mentioned before, since the effects of the Lee-Wick scalar are less than that of the normal scalar, it is still possible for the bounce to happen. Fig. 4 shows the equation of state of the system. We can see that the evolution of w is about the same as that in Fig. 1, since the change of the sign of d does not alter the result too much. Fig. 5 is the behavior of scale factor in this case while Fig. 6 gives the comparison of the energy densities of all components.
A change in the phase of the initial radiation field ve-  locity will not change the results (if we keep the other initial conditions fixed). On the other hand, if we flip the sign of the initial velocity of one of the two scalar fields, then the sign of the dominant contribution to the energy density as we approach the bounce will flip and this will prevent a bounce. If we use initial conditions containing two excited Fourier modes, then we obtain a bounce only if the signs of the initial field velocities are both chosen as in the first run whose results are shown here. Different phases for the scalar field velocities of the two modes destroys the possibility of obtaining a bounce.  Figures 7, 8 and 9 show the results for the equation of state parameter w, the Hubble parameter H and the contribution of the various components to the total ρ in the case of a simulation in which two Fourier modes are excited, with velocities of both scalar fields having opposite signs from those in the previous example. As is obvious, a Big Crunch singularity occurs. The scale factor of the universe driven by the system with parameters chosen as in Fig. 4. From the plot we see that the bounce happens at t = 0.  Fig. 4. The curves from up to down are: ρ φ 1 (black), ρ φ 2 (red), ρc (dark cyan) and ργ (blue), respectively. The variables are also normalized with the mass scale mrec = 10 −6 m P l .

VII. CONCLUSIONS AND DISCUSSION
In this paper we analyzed in detail the possibility of obtaining a cosmological bounce in a model which corresponds to the scalar field sector of the Lee-Wick theory coupled to relativistic radiation. It is known that the scalar field sector of the Lee-Wick theory in the absence of other fields can yield a cosmological bounce [1]. In fact, the universe will scale as non-relativistic matter with < w >≃ 0 both before and after the bounce. Thus, this model is a possible realization of the "matter bounce" scenario. However, this background is unstable to the introduction of radiation since in the contracting phase the growth of energy density in radiation will exceed that of matter and will lead to a Big Crunch singularity As has been shown in previous work [15], the introduction of a Lee-Wick partner to radiation does not prevent this instability. In this paper, we introduced an interaction between the radiation field and the scalar fields. The interaction could help drain energy from the radiation The one is of which the wavenumber k ≃ 0.01hM pc −1 with initial conditions γi ≃ −1.85 × 10 5 mrec,γi ≃ 7.35 × 10 6 m 2 rec , φ1i ≃ 1.015 × 10 5 mrec,φ1i ≃ 6.39 × 10 5 m 2 rec , φ2i ≃ 2.54 × 10 2 mrec,φ2i ≃ −4.96 × 10 3 m 2 rec , which if taken alone will give the bounce as has been shown in the previous example. The other is of which the wavenumber k ≃ 0.04hM pc −1 with initial conditions of the same initial values of the fields but the opposite signs of the scalar field velocity. From this plot we can see that the combination of the two Fourier mode will (generally) cause w blow up, thus preventing the bounce. This means that the bounce requires special fine-tuning of the initial phases for each Fourier mode. The coefficients c and d are chosen to be c = 10 −10 M 2 rec and d = −10 −10 M 2 rec .  Fig. 7. From the plot we see that there is a singularity at t = 0.
field to the Lee-Wick scalar, and thus could prevent the radiation from growing too fast to destroy the bounce. We analyzed the equations describing the evolution of the three matter fields (regular scalar field, its Lee-Wick partner and the radiation field) on a cosmological background both analytically and numerically. Our analytical analysis was perturbative and made use of the second order using Born approximation. The expansion parameter is set by the initial amplitude of the gauge field. We  Fig. 7. The curves from up to down are: ρ φ 1 (black), ρ φ 2 (red), ρc (dark cyan) and ργ (blue), respectively. The variables are also normalized with the mass scale mrec = 10 −6 m P l .
solved the equations of motion for each field at each order, and obtained their approximations in different cases. We compared their contributions to the total energy density, and derived necessary conditions for a bounce to happen. To support our analysis, we also performed numerical calculations.
Specifically, we investigated initial conditions in which one or two Fourier modes of the radiation field and the scalar field fluctuations are excited. We found special initial conditions which indeed lead to a non-singular bounce. Changing the sign of the initial scalar field velocity will destroy the bounce solution. In the presence of two Fourier modes, we found that a bounce requires identical initial phases for the two modes. For general initial conditions, we conjecture that the measure of such initial conditions which lead to a bounce is very small. We thus find that the addition of coupling terms between the scalar fields and radiation cannot save the Lee-Wick bounce background from the instability problem with respect to the addition of radiation (nor, for that matter, with respect to scalar field fluctuations). The instability problem with respect to anisotropic stress will be even worse.