text:
# Few Cycle Pulse Propagation

## I Introduction

The analysis of optical pulse propagation traditionally involves describing a pulse in terms of a complex field envelope, while neglecting the underlying rapid oscillations at its carrier frequency. The resulting “slowly varying envelope approximation” (SVEA) (see e.g. Shen (1984)), which reduces second order differential equations to first order, is valid when the envelope encompasses many cycles of the optical field and varies slowly. The alternative approach is to solve Maxwell’s equations numerically (see e.g. Brabec and Krausz (1997); Tarasishin et al. (2001)), which is more general but involves greater computational effort, and lacks the intuitive picture of a pulse “envelope”.

For example, optical parametric oscillators (OPOs) based on aperiodically-poled lithium niobate (APPLN) have generated 53 fs idler pulses at 3$`\mu `$m that are nearly transform limited, and contain only five optical cycles Beddard et al. (2000); laser pulses with less than three optical cycles have been generated in other contexts Baltuska et al. (1997). Under these circumstances, the validity of the slowly-varying envelope approximation is clearly open to question.

Brabec and Krausz Brabec and Krausz (1997) derived corrections to the SVEA, which they included in their “slowly evolving wave approximation” (SEWA). This enabled the few-cycle regime to be modelled with improved accuracy, and the SEWA has subsequently been applied in different situations, including ultrashort IR laser pulses in fused silica Ranka and Gaeta (1998); Tzortzakis et al. (2001), the filamentation of ultra-short laser pulses in air Akozbek et al. (2001), and even in microstructured optical fibres Gaeta (2002). Later, Porras Porras (1999) proposed a slightly different “slowly evolving envelope approximation” (SEEA) that included corrections for the transverse behavior of the field.

Here we use a field envelope approach to simulate the propagation of ultrashort pulses in a $`\chi ^{(2)}`$medium. The novelty is that we (a) derive a more general form than that of Brabec and Krausz, called the “generalised few-cycle envelope approximation”(GFEA); and (b) apply it to both optical (non-degenerate) parametric amplification (NPA), and the optical parametric oscillator (OPO). A shorter version of this paper has been published in Phys. Rev. A Kinsler and New (2003). The only previous multiple field application of this kind of result was for four wave mixing Gallagher-Faeder and Jonas (2000).

We compare the SEWA/SEEA equations to our own (section II), and explain the differences and subsequent adjustments to the necessary approximations. This theory enables us to rigorously study what combination of approximations affords the most efficient method for treating a given nonlinear interaction involving few-cycle pulses. Next (section III) we discuss the $`\chi ^{(2)}`$ nonlinearity and a scaling scheme designed to reveal the few-cycle effects. Then we compare the SVEA predictions to the few cycle GFEA theory using idealised situations (section IV) and more realistic OPO models (section V). Finally, we present our conclusions (section VI).

## II Theory

This section contains a summary of a complete rederivation Kinsler (2002) of a Brabec & Krausz style theory, which yields an evolution equation for an envelope description of pulse propagation in the few-cycle regime. Our result is more complicated than the SEWA equation Brabec and Krausz (1997), but reduces to it in the appropriate limits; it also explains the slight differences between their result and that of Porras Porras (1999).

Following Brabec-Krausz, we consider the case of small transverse inhomogeneities of the polarization, and so start with the three dimensional wave equation

$`\left(_z^2+_{}^2\right)E(\stackrel{}{r},t){\displaystyle \frac{1}{c^2}}_t^2{\displaystyle _{\mathrm{}}^t}𝑑t^{}ϵ(tt^{})E(\stackrel{}{r},t^{})`$
$`={\displaystyle \frac{4\pi }{c^2}}_t^2P_{nl}(\stackrel{}{r},t).`$ (1)
Here $`_{}^2`$ is the transverse Laplace operator, $`_\alpha `$ is shorthand notation for $`/\alpha `$, $`ϵ(t)=(2\pi )^1_{\mathrm{}}^{\mathrm{}}𝑑\omega \stackrel{~}{ϵ}(\omega )e^{ı\omega t}`$, $`\stackrel{~}{ϵ}(\omega )=1+4\pi \chi (\omega )`$, and $`\chi (\omega )`$ is the linear electric susceptibility. The electric field $`E`$ propagates along the $`z`$ direction. Both $`E`$ and the nonlinear polarization $`P_{nl}`$ are polarized parallel to the $`x`$ axis.

We can transform eqn. (1) into frequency space in order to expand $`\stackrel{~}{ϵ}(\omega )`$ in powers of $`\omega `$, thus enabling us to treat the material parameters as a power series which we can truncate to an approriate order. However for simplicity it is better to expand $`k`$ about a suitable $`\omega _0`$ instead. Using $`\stackrel{~}{ϵ}(\omega )=c^2k(\omega )^2/\omega ^2`$, it follows that

$`k(\omega )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\gamma _n\left(\omega \omega _0\right)^n}{n!}};`$
$`\gamma _n=_\omega ^nk(\omega )|_{\omega _0}=\beta _n+ı\alpha _n;\beta _n,\alpha _n.`$ (2)
We can now write the frequency space version of eqn. (1) as

$`\left(_z^2+_{}^2\right)E(\stackrel{}{r},t)+\left[{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ı^n\gamma _n\left(_t+ı\omega _0\right)^n}{n!}}\right]^2E(\stackrel{}{r},t)`$
$`={\displaystyle \frac{4\pi }{c^2}}_t^2P_{nl}(\stackrel{}{r},t).`$ (3)
We introduce an envelope and carrier form for the field in the usual way, using $`\stackrel{}{r}(\stackrel{}{r}_{},z)`$, so that

$`E(\stackrel{}{r},t)=`$ $`A(\stackrel{}{r}_{},z,t)e^{ı\mathrm{\Xi }}+A^{}(\stackrel{}{r}_{},z,t)e^{ı\mathrm{\Xi }},`$ (4)
and similarly $`P_{nl}(\stackrel{}{r},t)=B(\stackrel{}{r}_{},z,t;A)e^{ı\mathrm{\Xi }}+B^{}(\stackrel{}{r}_{},z,t;A)e^{ı\mathrm{\Xi }}`$. The symbol $`\mathrm{\Xi }=\beta _0z\omega _0t+\psi _0`$ is introduced as a convenient shorthand for the argument of the exponential. With these envelope-carrier substitutions, the equation of motion can be written as

$`\left(\left[ı\beta _0+_z\right]^2+_{}^2\right)A(\stackrel{}{r}_{},z,t)`$ (5)
$`+`$ $`\left[{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\gamma _n\omega _0^n}{n!}}\left({\displaystyle \frac{ı}{\omega _0}}_t\right)^n\right]^2A(\stackrel{}{r}_{},z,t)`$
$`=`$ $`{\displaystyle \frac{4\pi \omega _0^2}{c^2}}\left(1+{\displaystyle \frac{ı}{\omega _0}}_t\right)^2B(\stackrel{}{r}_{},z,t;A).`$
Eqn. (5) has the opposite sign on the RHS to Brabec & Krausz’s eqn. (2), but agreement is recovered later in eqn. (11).

As is usual, we introduce co-moving variables

$`\tau `$ $`=`$ $`\omega _0\left(t\beta _1z\right),_t\omega _0_\tau ,`$ (6)
$`\xi `$ $`=`$ $`\beta _0z,_z\beta _0_\xi \omega _0\beta _1_\tau ,`$ (7)
and eqn. (5) now becomes

$`\left\{\left(ı\beta _0+\beta _0_\xi \omega _0\beta _1_\tau \right)^2+_{}^2+\left[{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\gamma _n\omega _0^n}{n!}}\left(ı_\tau \right)^n\right]^2\right\}`$
$`\times A(\stackrel{}{r}_{},\xi ,\tau )+{\displaystyle \frac{4\pi \omega _0^2}{c^2}}(1+ı_\tau )^2B(\stackrel{}{r}_{},\xi ,\tau ;A)=0.`$ (8)
For convenience we also introduce the dimensionless ratio of phase and group velocities $`\sigma =\omega _0\beta _1/\beta _0=\left(\omega _0/\beta _0\right)/(1/\beta _1)=v_f/v_g`$, and use the fact that the refractive index at $`\omega _0`$ is $`n_0=c\beta _0/\omega _0`$. We also define a dispersion term $`\widehat{D}`$ in a similar way to Brabec-Krausz, but instead use a scaled (dimensionless) version $`\widehat{D}^{}=\left(\omega _0/\beta _0\right)\widehat{D}`$ in following equations so that

$`\widehat{D}^{}`$ $`=`$ $`{\displaystyle \frac{\omega _0}{\beta _0}}\left[ı\alpha _1\left(ı_\tau \right)+{\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\gamma _n\omega _0^{n1}}{n!}}\left(ı_\tau \right)^n\right].`$ (9)
Hence we get

$`0`$ $`=`$ $`\{(_\xi \sigma _\tau )+{\displaystyle \frac{1}{2ı}}(_\xi \sigma _\tau )^2+{\displaystyle \frac{1}{2ı\beta _0^2}}_{}^2`$ (10)
$`\left[ı\sigma \left(ı_\tau \right){\displaystyle \frac{ı\alpha _0}{\beta _0}}+ı\widehat{D}^{}\right]`$
$`+{\displaystyle \frac{ı}{2}}[ı\sigma \left(ı_\tau \right){\displaystyle \frac{ı\alpha _0}{\beta _0}}+ı\widehat{D}^{}]^2\}A(\stackrel{}{r}_{},\xi ,\tau )`$
$`+{\displaystyle \frac{2\pi }{ın_0^2}}\left(1+ı_\tau \right)^2B(\stackrel{}{r}_{},\xi ,\tau ;A).`$
This form can be rearranged without approximation to

$`_\xi A(\stackrel{}{r}_{},\xi ,\tau )`$ (11)
$`=`$ $`\left({\displaystyle \frac{\alpha _0}{\beta _0}}+ı\widehat{D}^{}\right)A(\stackrel{}{r}_{},\xi ,\tau )+{\displaystyle \frac{\left(ı/2\beta _0^2\right)_{}^2}{\left(1+ı\sigma _\tau \right)}}A(\stackrel{}{r}_{},\xi ,\tau )`$
$`+`$ $`{\displaystyle \frac{2ı\pi }{n_0^2}}{\displaystyle \frac{\left(1+ı_\tau \right)^2}{\left(1+ı\sigma _\tau \right)}}B(\stackrel{}{r}_{},\xi ,\tau ;A)+{\displaystyle \frac{T_R}{1+ı\sigma _\tau }},`$
where

$`T_R`$ $`=`$ $`\left[{\displaystyle \frac{ıq^2}{2}}_\xi ^2+{\displaystyle \frac{ı}{2}}\left({\displaystyle \frac{\alpha _0}{\beta _0}}ı\widehat{D}^{}\right)^2\right]A(\stackrel{}{r}_{},\xi ,\tau ).`$ (12)
Eqn. (11) is exact – it contains no more approximations than our starting point eqn. (1) except for the expansion of $`ϵ`$ in powers of $`\omega `$. We recover the full field $`E`$ from eqn. (4) by recombining $`A`$ and knowing the carrier. The partial derivatives ($`ı_\tau `$) in the denominators can, if necessary, be treated by Fourier transforming into the conjugate frequency space ($`\mathrm{\Omega }`$). Note that like $`\tau `$, $`\mathrm{\Omega }`$ is scaled relative to the carrier frequency.

If we set $`T_R=0`$, this gives us a generalised few cycle envelope (GFEA) equation, which contains the SVEA Shen (1984), SEWA Brabec and Krausz (1997), and SEEA Porras (1999) within it as special cases. Of course we cannot just set the $`T_R`$ term to zero without some justification, and this is discussed below.

The $`\left(2ı\pi /n_0^2\right)𝒦B`$ polarization term from eqn. (11) has prefactors which depend on the time-derivative of the polarization, and these new terms are what add the effect of finite pulse lengths to the pulse evolution. Note that we can write this polarization term in different forms:

$`𝒦`$ $`=`$ $`{\displaystyle \frac{\left(1+ı_\tau \right)^2}{\left(1+ı\sigma _\tau \right)}}`$ (13)
$`=`$ $`\left(1+ı\sigma _\tau \right)\left[1+\left(1\sigma \right){\displaystyle \frac{2ı_\tau +\left(1+\sigma \right)_\tau ^2}{\left(1+ı\sigma _\tau \right)^2}}\right]`$
$`=`$ $`\left(1+ı_\tau \right)\left[1+{\displaystyle \frac{ı\left(1\sigma \right)_\tau }{\left(1+ı\sigma _\tau \right)^2}}\right].`$
With $`\sigma =1`$, these reduce to the $`𝒦=1+ı_\tau `$ SEWA Brabec and Krausz (1997) form. Similarly, to first order in $`\left(\sigma 1\right)`$, one can get the $`𝒦=1+ı\sigma _\tau `$ form analogous to the SEEA corrections Porras (1999), although a strict expansion to first order in $`_\tau `$ Kinsler (2002) gives a prefactor of $`\left[1+ı\left(2\sigma \right)_\tau \right]`$. Finally, for a SVEA theory, $`𝒦=1`$, since the $`_\tau `$ terms are assumed to be negligible.

The $`T_R`$ term is negligible if the following conditions hold:

terms in $`_\tau `$ can be neglected if

$`\left|\left({\displaystyle \frac{\omega _0^m\gamma _m^{}}{\beta _0m!}}\mathrm{\Omega }^m\right)\stackrel{~}{A}(\stackrel{}{r}_{},\xi ,\mathrm{\Omega })\right|`$ $`\left|\stackrel{~}{A}(\stackrel{}{r}_{},\xi ,\mathrm{\Omega })\right|`$ (14)
terms in $`_\xi ^2`$ can be neglected if

$`\left|_\xi \stackrel{~}{A}(\stackrel{}{r}_{},\xi ,\mathrm{\Omega })\right|\left|\stackrel{~}{A}(\stackrel{}{r}_{},\xi ,\mathrm{\Omega })\right|,`$ (15)
and eqn. (15) only holds if, in addition,

terms in $`_{}^2`$ can be neglected if

$`\left(1+\sigma \mathrm{\Omega }\right)\beta _0^2w_0^21,`$ (16)
is “weak” if

$`{\displaystyle \frac{n_0^2}{2\pi }}{\displaystyle \frac{\left(1+\sigma \mathrm{\Omega }\right)}{\left(1+\mathrm{\Omega }\right)^2}}{\displaystyle \frac{\left|\stackrel{~}{B}(\stackrel{}{r}_{},\xi ,\mathrm{\Omega };A)\right|}{\left|\stackrel{~}{A}(\stackrel{}{r}_{},\xi ,\mathrm{\Omega })\right|}}.`$ (17)
We use $`\mathrm{\Omega }`$ instead of $`ı_\tau `$ for these conditions because conditions on the frequency components of the various terms are better defined than those for time derivatives.

These conditions are the same as those required for the SEWA and SEEA theories, with the SVEA conditions being a special case gained by setting $`\left|\mathrm{\Omega }\right|1`$ for the diffraction and nonlinearity conditions – implying that modulations in the envelope are so slow compared to the carrier frequency that they can be neglected. Note that backwardly propagating behaviour has not been explicitly excluded, but since it would appear as a modulation on the envelope $`A`$, it would be approximated away as part of the evolution condition (eqn. (15)).

Note that the exact solution of eqn. (11) makes no reference to a particular choice of carrier phase $`\psi _0`$. This implies that once a solution for the propagation of a particular envelope has been obtained, the problem has in fact been solved for a set of pulses (and initial conditions) based on different carrier phases – where that set is determined by the initial envelope and some arbitrary choices of carrier phase $`\psi _1[0,2\pi )`$. The final state is then given by the chosen $`\psi _1`$ combined with the final form of the envelope.

## III The $`\chi ^{(2)}`$ nonlinear system

When modelling $`\chi ^{(2)}`$ nonlinear systems we split the optical field into two or three parts, depending on whether a degenerate or non-degenerate system is being treated. For example, a parametric amplifier would have pump, signal, and idler field components. We then define an envelope $`A_\alpha `$ and carrier $`e^{ı\mathrm{\Xi }_\alpha }`$, $`\mathrm{\Xi }_\alpha =\beta _{\alpha ,0}z\omega _{\alpha ,0}t+\psi _{\alpha ,0}`$ for each field component, and use a separate propagation equation for each. The total field is then the sum of these different components:

$`E`$ $`=`$ $`E_p+E_s+E_i`$ (18)
$`=`$ $`A_pe^{ı\mathrm{\Xi }_p}+A_p^{}e^{ı\mathrm{\Xi }_p}+A_se^{ı\mathrm{\Xi }_s}+A_s^{}e^{ı\mathrm{\Xi }_s}`$
$`+A_ie^{ı\mathrm{\Xi }_i}+A_i^{}e^{ı\mathrm{\Xi }_i}`$
Because the wave equation eqn. (1) is linear in the electric field, we can use eqn. (18) in the theory of section II, and get a GFEA-like equation which is rather like the sum of three separate copies of eqn. (11), but with a single polarization term. In a $`\chi ^{(2)}`$ medium this polarization term is proportional to the square of the field, and inspection of eqn. (18) shows that this can be expanded into a sum of many terms. We can then split the multiple-field GFEA equation into three separate GFEA-like equations, one for each envelope. In doing this we need to make sure to assign the pieces of the polarization term appropriately to a suitable envelope equation, whilst discarding those that are not resonant with a field carrier.

Our chosen nonlinear crystal is congruent LiNbO3, for which we calculate refractive index and dispersion data from the Sellmeier equations of Jundt Jundt (1997). We model the nonlinear polarization using the square of the total electric field, retaining the parts resonant with our field carriers in the normal way. Our OPO simulations (see section V) assumed a pump frequency of 357.1 THz, with nominal signal and idler carrier frequencies of 257.5 THz and 99.6 THz respectively (wavelengths $`0.84000\mu `$m, $`1.16500\mu `$m, $`3.0110\mu `$m). This means the pump pulse will travel through LiNbO3 more slowly than the signal and it therefore needs to be injected into the crystal ahead of it (see fig. 1). When the two overlap, an idler is generated by the nonlinear interaction, and the three pulses then continue to interact with each other as they propagate through the crystal. Note that our ideal non-degenerate parametric amplifier simulations (see section IV) use the same field frequencies, but idealise the crystal parameters by setting the group velocities and dispersions to zero.

### III.1 System Scalings

In a typical experiment, the crystal length would be fixed, as would any properties defined by its design, such as periodic poling. If we were to investigate this case for a range of pulse durations, there would naturally be differences between the results, even within the SVEA. For example, the relative pulse broadening caused by travelling through a $`1000\mu `$m crystal is greater for a 12fs pump pulse than for a 48fs one. Similarly, a fixed timing offset for injection would have different effects; and a fixed pump pulse power would generate different strengths of nonlinear interaction for different pulse lengths. All these effects would confuse any attempt at a systematic comparison of the few cycle effects in the models we consider.

Therefore, in order to isolate specific few-cycle effects, we must scale the pump pulse FWHM $`𝒯_\text{p}`$, crystal length $`L_\text{c}`$, pump pulse energy $`𝒲`$, and pump timing offset $`t_\mathrm{\Delta }`$ in such a way as to ensure the effects of group velocity, pump timing, and nonlinearity occur in the same proportions to one another over the range of pulse lengths.

We can work out an appropriate scaling by examining a simple version of the propagation equation (eqn.11), where we write the group velocity prefactors as $`B_1`$, the second order dispersion prefactors as $`B_2`$, and the polarization terms as $`CA^2`$. To assist us with the scaling process we also write $`\xi =r^f\xi ^{}`$, $`\tau =r^g\tau ^{}`$, and $`A=r^hA^{}`$, where $`r`$ is the scaling factor. Our simple propagation equation is therefore

$`r^{h+f}_\xi ^{}A^{}`$ $`=`$ $`r^{h+g}B_1_\tau ^{}A^{}+r^{h+2g}B_2_\tau ^{}^2A^{}+r^{2h}CA^2.`$ (19)
We can easily match the LHS term with the polarization term by setting $`f=h`$; but then we must choose either $`h=g`$ to match group velocities, or $`f=2g`$ to match the second order dispersion – we cannot match both. For our chosen OPO situation (see section V), it is best to match the group velocity terms, which control how long the pump and signal pulses overlap – in general, the dispersion has a much smaller effect.

We take our reference situation to be a $`20`$nJ $`24`$fs FWHM pump pulse propagating through a $`500\mu `$m crystal, with a pump timing offset of $`48`$fs. For the chosen parameter scaling ($`f=g=h`$)

$`{\displaystyle \frac{𝒯_\text{p}}{24\text{fs}}}={\displaystyle \frac{L_\text{c}}{500\mu \text{m}}}={\displaystyle \frac{t_\mathrm{\Delta }}{48\text{fs}}}={\displaystyle \frac{20\text{nJ}}{𝒲}}.`$ (20)
This removes the gross effects caused by reducing the pulse duration, so that if there are no few cycle effects, and second order dispersion is negligible, each pulse length should give identical results.

First, the crystal length $`L_\text{c}`$ is scaled in direct proportion to the pulse length $`𝒯_\text{p}`$. This means that the relative group velocities are scaled so that the pump and signal pulses spend the same relative time overlapping, leading to an equivalent nonlinear effect.

Second, and related to the previous point, the pump timing offset $`t_\mathrm{\Delta }`$ is reduced in direct proportion to the crystal length (and hence pulse length). This means the pump-signal pulse cross-over occurs at the same relative point, so the effect of the pump-signal overlap is similar. If a 48fs pump pulse in a $`1000\mu `$m crystal arrives $`96`$ns before the signal, the crosover occurs at about the $`500\mu `$m half-way mark – just as for a 12fs pulse in a $`250\mu `$m crystal, where a $`t_\mathrm{\Delta }=24fs`$ gives a crossover at the $`125\mu `$m half-way mark.

Third, a shorter crystal reduces the effect of the nonlinearity, so the pump pulse energy is scaled up as the inverse of the reduction in crystal length – a half length crystal needs pulse energies $`𝒲`$ twice as large – a 10nJ 48fs pulse in its $`1000\mu `$m crystal with see the same nonlinear effect as a 20nJ 24fs pulse in its $`500\mu `$m crystal. Note that the pulse energy depends on the area of the pulse, and so depends on its length, so that $`𝒲\left|A\right|`$.

We could choose to make the scaling perfect, by also scaling the crystal parameters. If we scale the crystal dispersion with $`B_2=r^gB_2^{}`$, the relative amount of pulse spreading changes to become the same for each simulation – e.g. if the 48fs pulse widens by 10% in a $`1000\mu `$m crystal, the 12fs will also widen by 10% in its $`250\mu `$m crystal. We did a set of SVEA simulations on this basis, and as expected saw identical pulse profiles regardless of the chosen pulse length. However, we chose not to use this perfect scheme for the bulk of our OPO simulations because it is far from being experimentally practical.

## IV Ideal Parametric Interactions

A parametric amplifier is a single-pass device: pump and signal pulses are injected into one end of the nonlinear crystal, they interact within it, then exit at the far end. However, because real nonlinear crystals (such as LiNbO3) tend to have significant dispersion, very short pulses quickly spread out, making them difficult to create, and reducing the few-cycle effects we aim to study.

In order to demonstrate clearly the nature of few-cycle effects in $`\chi ^{(2)}`$ materials, in this section we investigate an ideal case by setting the dispersion to zero, and make the group velocity the same at all frequencies. This means that $`\sigma =1`$, so the “few-cycle” polarization prefactor $`𝒦`$ is identical for both the SEWA and GFEA theories. Note that it is difficult to do no-dispersion simulations over long times because pulse self-steepening causes both the numerical integration and the theoretical approximations to break down.

We inject Gaussian pump and signal pulses at exactly the same time (i.e. $`t_\mathrm{\Delta }=0`$), with the same width. They then travel down the crystal with maximum overlap, interacting all the way. Other parameters are fixed by the scaling rules from section III.1. Further, when graphing results for the figures, we scale the times for each pulse length to the 6fs case (e.g. for a 24fs pulse, “$`\tau =10`$” corresponds to $`40`$fs), and scale the pulse intensities in proportion to their initial intensities. This means that graphs of the initial conditions for a range of pulse lengths would be identical.

Finally, note that in these ideal results, the nonlinear interaction is “strong”, with significant transfer of energy between the fields.

### IV.1 Second Harmonic Generation (SHG)

In second harmonic generation the nonlinearity causes a field at one frequency (the fundamental) to generate one at twice that frequency (the second harmonic). In our 48fs simulations, a 10nJ pulse is injected in at the fundamental frequency (see fig. 2).

In Figs. 4, 5 we can see comparisons of the pulses resulting from a single pass through the crystal. Although the plotted (and scaled) $`\left|A\right|^2`$ of the pulses are rather similar; the GFEA theory gives noticeably different results for the shortest pulses, unlike the SVEA theory. Note that the vertical scales on the two figures (not shown) are rather different.

However, this apparent similarity in $`\left|A\right|^2`$ is a little misleading because the underlying phase profiles of the envelopes are very different, as can be seen in (see fig. 6). More features, primarily extra oscillations, do develop at shorter pulse durations, but generally the comparisons remain similar to those shown.

### IV.2 Degenerate Parametric Amplification (DPA)

A degenerate parametric amplifier uses energy from a pump pulse to amplify that of a signal pulse at half the frequency ($`\omega _p=2\omega _s`$, see fig. 7). For our 24fs reference simulations, a 20nJ pulse is injected at the pump frequency, and a signal pulse of 10pJ is amplified; the nonlinear interaction is very strong.

In Figs. 9, 10 we can see comparisons of the pulses resulting from a single pass through the crystal. Note that the intensity profiles of the pulses are very similar, but the GFEA theory gives different results to the SVEA for the shortest pulses.

However, this apparent similarity in $`\left|A\right|^2`$ is a little misleading because the underlying phase profiles of the envelopes are very different, as can be seen in fig. 11. The change in phase of the pump pulse in both the SVEA and GFEA case is due to the presence of nodes that develop to either side of its centre. More features, primarily extra oscillations, do develop at shorter pulse durations, but generally the comparisons remain similar.

### IV.3 Degenerate Parametric De-amplification (DPD)

This is the same situation as for degenerate parametric amplification, but we have changed the relative phase between the pump pulse and signal so that the input signal pulse is de-amplified (see fig. 12). We are interested in this case because the de-amplification is very sensitive to the phase relationship between the pulses, and the few-cycles terms modify the phases of the pulses. Our SVEA model of the situation behaves in the expected way, with the signal pulse decaying away to zero as it propagates through the crystal and transfers its energy to the pump.

However the new terms in the GFEA propagation equation adjust the phases of both the pump and signal fields as they propagate, moving them away from an exact match to the de-amplification criteria. As the mismatch increases, the de-amplification process is reversed, and instead the usual amplification takes over. Fig. 14 shows the output signal pulse profiles on exit from the crystal, note the enormous difference in the few-cycle dependence of the SVEA and GFEA theories – both in character and amplitude. The SVEA profile is just the residual input signal pulse which has not yet been fully “de-amplified”; this residue can also be seen as a shoulder on the wings of the GFEA profiles if plotted on a $`\mathrm{log}`$ scale. The large and very visible double peak in the GFEA graph is a result of the finite pulse length on the nonlinear polarizations term – as can be seen in eqn. (13), this alters the phase profile of the nonlinear effect, and hence the pulses, which means the signal no longer satisfies the exact criteria for de-amplification. For our chosen input pulse powers, this effect persists for initial phase mismatches up to $`0.01`$ radians, after which the amplification of the signal pulse has largely swamped the anisotropy induced by the few-cycle effects. Of course the effect is most visible when the interaction has been strong enough for the input component of the signal pulse to (almost) completely disappear.

In fig. 16 we see how the behaviour changes both with pulse length and initial phase match – we see perfect de-amplification in the SVEA model, but the steeply sloping $`\varphi _{\mathrm{signal}}=\pi /4`$ line clearly demonstrates the effects of a finite pulse length. Of course as the $`\varphi _{\mathrm{signal}}`$ moves further from $`\pi /4`$, the SVEA model no longer undergoes perfect de-amplification, and so the difference between the two models becomes small.

### IV.4 Non-degenerate Parametric Amplification (NPA)

We consider first a non-degenerate parametric amplifier with pump, signal, and idler frequencies such that $`\omega _p\omega _s+\omega _i`$, and $`\omega _s\omega _i`$. In the 24fs reference case, the initial pump energy is 20nJ and the initial signal energy is 10pJ, with a negligible (but finite) idler. For other pulse durations, the energies were scaled according to eqn. (20). Fig. 17 shows how, according to the GFEA, the idler pulse intensity profiles $`\left|A_i\right|^2`$ generated in a single pass of the crystal vary with pulse duration. The profiles show little variation with pulse duration except for the shortest pulses ($`\tau 20`$), where distortion is evident; the signal and pump profiles show deviations of a comparable magnitude.

In fig. 19, 20, 21 we can see equivalent comparisons for the signal and pump intensity profiles.

More dramatic effects appear in the phase profiles: in fig. 22, the phases of the pulse envelopes at pulse durations of 18fs and 96fs are shown with the phase distortions due to the finite pulse lengths (see eqn. (13)). As the pulse duration shortens, the principal effect is to increase the magnitude of the phase distortion, leaving the shape of each profile largely unchanged; however more complex phase oscillations develop for the shortest pulses. At 96fs, the profiles show a smaller distortion, and are tending towards the long-pulse SVEA limit. In this limit, the profiles are essentially flat, although the pump field develops nodes which give rise to a step-like change in the phase.

### IV.5 Non-degenerate Parametric De-amplification (NPD)

As a variant on the case just treated, signal and idler pulses with equal numbers of photons were injected, and the relative phases of the pulses set to ensure that the signal and idler experience initial de-amplification (see fig. 23). Since the subsequent evolution is sensitive to phase changes, and the finite pulse length terms in the GFEA affect the phases, this is an interesting situation to examine. In the SVEA, the signal and idler decay away towards zero as the pulses propagate, so the SVEA output signal is just some residual part of the input. The GFEA evolution is different, as can be seen from eqn. (13) – the finite pulse lengths alter the phase profile of the nonlinearity, and hence change the evolution of the pulses. During an initial period of de-amplification, the pulses undergo a gradual phase distortion. Then, as the discrepancy increases, amplification takes over. In a comparison of SVEA and GFEA models, the effect caused by the phase distortion is more visible when the interaction is been strong enough for the input component of the signal pulse to be strongly depleted, and also is much stronger for shorter pulses

The GFEA signal pulse profiles on exit from the crystal as a function of pulse duration are presented in fig. 25. Note that the SVEA prediction corresponds to the long-pulse limit of the GFEA figure, but those limiting features are too small to be seen.

The double peak is a result of the finite pulse length on the nonlinear polarization – the “phase twist” caused by the few cycle terms is odd, causing a node at $`\tau =0`$ as the field is amplified – at exactly $`\tau =0`$, the field continues to be de-amplified. Note that at 18fs, for example, the output signal pulse amplitude is larger than the input value.

The SVEA profile in fig. 24 is just the residual input signal pulse which has not yet been fully “de-amplified”; this residue can also be seen as a shoulder on the wings of the GFEA profiles if plotted on a $`\mathrm{log}`$ scale. Fig. 25 shows the output GFEA signal pulse profiles on exit from the crystal, note the enormous difference in the few-cycle dependence of the SVEA and GFEA theories – both in character and amplitude.

In fig. 26 we see how the SVEA behaviour changes both with pulse length and initial phase – here we see perfect de-amplification.

The GFEA output pulse energies are displayed in fig. 27, which shows how the behaviour changes both with pulse length and initial phase. The data for $`\varphi _{\mathrm{signal}}=\pi /2`$ demonstrates the effects of exact initial conditions and finite pulse length; maximum de-amplification occurs in the long-pulse (SVEA) limit. If we instead start with a signal phase slightly different from $`\pi /2`$, e.g. $`0.51\pi `$, the de-amplification is less efficient and will eventually be overtaken by the amplification, even for the SVEA model. Consequently, comparisons for imperfect initial phases are dependent on the length of the crystal. However, since we use a scaling procedure, the results still behave in a systematic way, even if they are not completely generic.

Of course, changing other initial conditions can also disturb the de-amplification: e.g. different numbers of signal and idler photons. Although both signal and idler will initially be de-amplified, as they approach zero photon number, one field will “overshoot” the zero and be inverted. This alters the phase relationships, and so again amplification takes over. As an example, simulations based on our 18fs pulses suggested that photon number mismatches of about one percent would not noticeably disrupt the appearance of either fig. 25 or 27.

These comparisons testing the sensitivity of these results to mismatches in the initial conditions are of course dependent on our chosen system parameters. For example, either more intense pulses or longer crystals would diminish the difference between SVEA and GFEA predictions; conversely, weaker pulses or shorter crystals would enhance them.

## V Optical Parametric Oscillation (OPO)

We move on from optical parametric amplification to a synchronously-pumped Optical Parametric Oscillator (OPO). As shown in fig. 28, we considered the case of a LiNbO3 crystal in an optical cavity with mirrors that reflect the signal wavelength only. The oscillator is driven by a train of gaussian pump pulses whose periodicity closely matches the natural period of the cavity, and which amplify and then sustain the signal pulse confined within it. The cavity length can be “tuned” about exact synchronisation. The idler pulse, generated when the signal pulse interacts with each new pump pulse, is transmitted through the output mirror with the pump, while the signal is strongly reflected. For a given set of parameter values, we modelled the development of the signal pulse over many cavity transits until it reached a steady state. Typically, we found that the signal stabilised in several hundred transits although, in a few cases, no equilibrium was achieved and the system oscillated indefinitely. We studied the evolution for a range of phase-matching conditions ($`\mathrm{\Delta }k=0`$ to $`24\times 10^3/\mu `$m, and cavity length tunings ($`16`$fs to $`+16`$fs). Except where stated, the figures contain data for the perfectly phase matched and sychronised case.

The steady state was determined by propagating the pulses through many passes of the oscillator until the modulus-squared of the pulse envelopes had stabilized.

### V.1 Dispersion Scaled OPO

As discussed in section III.1, it is possible to make the SVEA propagation equation completely scale invariant if we modify the crystal dispersion parameters as well as the pulse length, energy, and crystal length. Although this is not an experimentally achievable goal, it is instructive to look at what happens to the OPO output in this instance. To make automatic processing of the results more convenient, we adjusted the crystal lengths to be exact powers of two, which differed slightly from those derived from the standard reference length: i.e. the 24fs pulse case related to a crystal length of 512$`\mu `$m, not 500$`\mu `$m. As expected from the scaling, and shown on figs. 29, 30, the scaled SVEA results were independent of pulse length.

In contrast the GFEA pulse profiles have a more complex behaviour, as shown on figs. 31, 32. For the chosen parameters, they have a profile which becomes increasingly distorted at shorter pulses; but there is also an abrupt transition to what is an SVEA-like pulse profile at between 36 and 48fs.

On fig. 33 we can see this transition and the expected gradually convergence of GFEA and SVEA; the maximum differences plotted are about 10% of the pulse heights, this message is repeated on fig. 34. If the convergence of GFEA to SVEA seems slower than inutitively expected, this could be partly explained by noting that the effect of the extra $`_\tau `$ term caused by the non slowly-varying nature of the pulses has less distance over which to accumulate, so a similar discrepancy between SVEA and SEWA at e.g. 24fs and 12fs means that the 12fs case accumulates additional evolution at twice the rate as at 24fs – as indeed is to be expected, since broadly speaking, the time derivative of a pulse envelope will double as its pulse length is halved.

Note: the slight disruption to the wings of the 92fs results are because the time width allowed for by the simulation is smaller when compared to the pulse width for this particular case.

### V.2 Scaled OPO

The complex nature of the dynamics, which arises from repetitive cycling of the signal pulse in the presence of many interacting processes makes the isolation and analysis of few-cycle effects within the different models quite complicated. Fig. 35 shows intensity profiles for the pump, signal, and idler (bottom to top in each frame) for the SVEA (dashed line) and GFEA (solid line) for four different pulse durations. Simulations with pump pulse durations of 48fs and over never reached a steady state, and so are not included here. The first point to note in fig. 35 is that the SVEA results are not identical in all frames, even though the scaling procedure in section III.1 is designed to make them, as far as possible, independent of pulse duration. The reason is that, as noted in section III.1, the dispersion scales in a different way to the group time delay, and so is not correctly compensated by eqn. (20).

A second rather surprising feature is that we might expect the GFEA results to tend to the SVEA as pulse length increases, but this is not evident from the graphs. The explanation for this is that the steady state of the OPO can change suddenly as the parameters are varied. This property is highlighted in fig. 31, which shows the GFEA signal pulse profile for pulse durations from 6fs to 192fs; the sudden adjustment of the GFEA when moving from 36fs to 48fs takes it close to the SVEA, and the difference between the two gradually disappears as the pulse duration is increased further (see Fig. 34 ). Note that the scaling procedure used for fig. 31 is an extension of eqn. (20) in that the dispersion term is also scaled, making the SVEA results completely independent of pulse duration.

Spectral profiles corresponding to the temporal profiles of fig. 35 and 36 are shown in fig. 37. The spectral shape for each field is similar across all pulse durations, with a pulse of double the (time) width naturally having half the bandwidth. Notice that the pump and signal spectra in the 6fs frame are close to overlapping, which indicates that the separation of the total EM field into distinct pump, signal, and idler components is becoming a questionable assumption.

Inclusion of the carrier wave in the results raises some quite subtle issues that need careful consideration. It must be stressed again that the carrier drops out of the analysis leading to eqn. (11). The envelope description is therefore complete, although the phases of two of the three envelope functions can be changed by arbitrary constants without any effect on the computations apart from an appropriate adjustment in the phase of the third envelope. For instance, if the phases of the pump and signal envelopes are changed by $`\mathrm{\Delta }\varphi _p`$ and $`\mathrm{\Delta }\varphi _s`$, the phase of the idler envelope is changed by $`\mathrm{\Delta }\varphi _i=\mathrm{\Delta }\varphi _p\mathrm{\Delta }\varphi _s`$. Adjustments of this kind show up in the results only if graphs of the complete electric field profiles, including the carrier waves, are displayed, as in fig. 38. If the simulations in that figure were re-run with differing envelope phases, this would be reflected in temporal displacements of the carrier-like oscillations beneath the envelopes.

A further interesting feature is that, while the moduli of the pulse enveleopes may have stabilised in a simulation, the envelope phases can (and usually do) change from pass to pass; this process continues indefinitely, so a movie made up of frames from successive transits would show the pump, signal, and idler electric field oscillations drifting across underneath the respective steady envelope profiles.

Note that although it is common to talk of “shifts in carrier phase”, in the theory the carrier phase is specified by the initial conditions, and does not change. The effects usually ascribed to “shifts in carrier phase” instead manifest themselves as a shift in the phase of the complex envelope.

The different models discussed in this paper give significantly different results for the pass-to-pass phase drift. Fig. 39 shows the phase change for the signal pulses as a function of pulse length for the SVEA, SEWA, and GFEA; note that the SEWA and GFEA results are similar to each other, while the (less accurate) SVEA exhibits a very different dependence.

The reference point used in calculating the phase drift is at the maximum amplitude of the envelope of the signal pulse, which is in fact not necessarily at the point of maximum electric field. This is a good choice for our purposes because it does not move between passes once a steady state is established. Although these phase drifts are quite small, discrepancies between the SVEA and GFEA will quickly accumulate.

Analysis shows that the phase of the envelope at this reference changes by a fixed amount from one pass to the next. Since this phase drift tends to be small, a slow evolution of the E field profile can be seen between passes in any given simulation. Further, the accumulated difference between the less accurate SVEA prediction and the GFEA can lead to large differences in the predicted electric field profiles, even if the envelope profiles happened to be similar. Note that the idler pulse envelope phase drift is in the opposite direction to the signal drift, since the pump pulse arrives with the same envelope phase at the beginning of each pass.

The details of fig. 39 are less important than its message. Getting an output of either an identical idler pulse from pass to pass or at least a well understood pulse shape progression can be important, and this will be achieved using accurate models or by experiment based on accurate models. We see that even for relatively long pulses the SVEA theory predicts a different phase drift to that of the SEWA/GFEA theory; and that the trend is different even when using a sequence of parameters designed to minimise the differences between simulations of different pulse lengths.

### V.3 Synchronisation and Phase Mismatch

The phase drift varies as a function of the phase mismatch and the cavity synchronisation: the data is shown in Figs. 40, 41. The glitches in the data are caused by occasional difficulties in evaluating the phase of the envelopes and the wide variations that can occur for near-zero fields – not all the data points represent good steady states.

### V.4 Fixed Length Crystal

Here we present results for a more realistic case of a fixed length crystal ($`1000\mu `$m) and a range of pump pulse durations ($`6`$$`96`$fs); with all parameters being kept constant, notably a fixed pulse energy. Since the crystal is relatively long, the high dispersion comes into play, ensuring that even for rather short pump pulses (e.g. $`6`$fs), the signal and idler become relatively broad ($`200`$fs) – thus the role of few-cycle effects should be relatively small. However, inspection of the envelope profiles for the SVEA and GFEA models still show small but noticeable differences, for both perfect and imperfect phase matching, even for pump pulses as long as $`96`$fs.

In fig. 42 we see how the inter-pass phase difference for the signal (and hence idler) for a fixed 1000$`\mu `$m crystal with decreasing pump pulse widths and adjusted time offsets, but fixed pulse energy. Difference between the SVEA predicted phase drift and the GFEA are noteable for perfect phase matching ($`\mathrm{\Delta }k=0`$), but when this is no longer exact, the two predictions rapidly become similar, as now both a long crystal and imperfect phase matching both act against the generation of clear few cycle effects.

## VI Conclusions

We have presented a new and more complete derivation of how the envelopes of extremely short optical pulses evolve in nonlinear interactions. We have compared the results of our new (GFEA) model to those of the traditional slowly varying envelope approximation (SVEA) using a scaling procedure to distinguish specific few-cycle effects from other phenomena caused by changing pulse duration. It should be noted that the SVEA becomes inadequate whenever the envelope changes rapidly within a few carrier periods. Strictly speaking, a few-cycle pulse is not required, because a steep edge within a longer pulse also fulfils the conditions.

The effect of the extra “few-cycle” terms in the GFEA evolution equation is to add a phase distortion to the nonlinear polarization term, which then imposes itself on the pulse envelopes. This is demonstrated by our single-pass optical parametric amplifier NPA model where, whilst the SVEA model is insensitive to pulse length, the GFEA theory shows clear changes as the pulses get shorter and contain fewer optical cycles.

Further, when we studied the highly sensitive de-amplification case (i.e. NPD), we saw dramatic differences between the SVEA and GFEA simulations even outside the few-cycle regime. These arose from the phase distorting effects of the few-cycle terms in the theory disrupting the exact phase relationships needed for de-amplification. While the absolute size of these differences do depend on the chosen parameters of crystal length, pulse energy, and so on, they will always get dramatically larger for shorter pulses,

On the other hand, the repetitive cycling nature of the optical parametric oscillator (OPO) produces more complicated and subtle dynamics; small changes in parameter values can, for instance, cause sudden changes in the steady state fields. It is therefore no surprise that comparison of the results predicted by the different models is less straightforward in the OPO case. The new model certainly produces differences in the pulse envelopes as well as the phases, although the way in which the GFEA tends to the SVEA in the long-pulse limit has some interesting features. The two models also predict different results for the pass-to-pass phase drift of OPO pulses, and this implies significant differences in the electric field structures. In both cases, the carrier wave moves under the envelope from one transit to the next, but by different amounts.

It can be useful to regard these “few-cycle” effects as adding a phase twist to the envelope evolution. This then shows itself most clearly in the pass-to-pass phase drift in the signal and idler phases (see Fig. 39). Finally, although these effects are usually described as “few-cycle effects”, they are more accurately described as “finite pulse length” effects. Although for the parameter ranges needed to describe typical nonlinear crystals, it is only in a few-cycle regime where the effects are easily visible in the pulse intensity profiles, these twisting effects on the phase structure of the pulse do occur, and could be seen in a many cycle regime.


format: text

text:
# 1 Introduction

## 1 Introduction

Dark energy or missing energy is one of the contemporary issues the physics community is more interested in due, mainly, to a relatively recent (revolutionary) discovery that our present universe is in a stage of accelerated expansion. This missing component of the material content of the universe is the responsible for the current stage of accelerated expansion and accounts for 2/3 of the total energy content of the universe, determining its destiny. This is a new form of energy with repulsive gravity and possible implications for quantum theory and supersymmetry breaking. Among others, the new cosmology is characterized by the following features: 1) Flat, critical density accelerating universe, 2) early period of inflation, 3) density inhomogeneities produced from quantum fluctuations during inflation, 4) composition of 2/3 dark energy, 1/3 dark matter and 1/200 bright stars, 5) matter content: $`(29\pm 4)\%`$ cold dark matter, $`(4\pm 1)\%`$ baryons and $`0.3\%`$ neutrinos. Besides, there is evidence (recent observation of SN1997ff at redshift $`z=1.7`$) that the present stage of accelerated expansion was preceded by an early period of decelerated expansion. Any theoretical framework intended to describe our present universe should be compatible with these observational evidences.

A self-interacting, slowly varying scalar field, most often called quintessence, has been meant to account for the dark energy component. In the simplest models of quintessence, minimally coupled scalar fields are considered. However, extended quintessence models have also been considered. In these models the scalar field that accounts for the dark energy component is non-minimally coupled to gravity and, besides, is a part of the gravity sector. Brans-Dicke theory is a prototype for this kind. Motivated by the possibility to have superaccelerated expansion and, correspondingly, superquintessence solutions, this type of theories is now being considered (see, for instance, Ref. and references therein). In particular, in Ref., the author considered different non-minimally coupled theories (among them the Brans-Dicke theory) and showed that, although superquintessence may arise in these models, its effect is vanishingly small, i. e., is far beyong the present observational possibilities. In that paper the author numerically integrated the field equations by using a computed code. Instead, it would be of interest to have a parametric class of exact solutions to treat this subject.

The aim of the present paper is, precisely, to derive a new parametric class of exact cosmological solutions to Brans-Dicke gravity theory with a self-interacting scalar field, by assuming a linear relationship between the Hubble expansion parameter and the time derivative of the scalar field. This type of relationship is a particular case of a more general one that has been used in Ref.  to find general solutions in $`N+1`$-dimensional theories of gravity with a self-interacting (minimally coupled) scalar field and also to derive general solutions when we have a self-interacting scalar field with exponential potential plus a free scalar field.111A method where a linear relationship between the field variables and/or their derivatives is assumed has also been used in  to derive 4d Poincare invariant solutions in thick brane contexts. Also, it has been already used by us to derive exact cosmological solutions to Einstein’s gravity with two fluids: a barotropic perfect fluid of ordinary matter, together with a self-interacting scalar field fluid accounting for the dark energy in the universe. The assumed relationship between the Hubble parameter and the time derivative of the scalar field is suggested by an implicit symmetry of the field equations.

We point out that it is not neccessary to make any a priori assumptions about the functional form of the self-interaction potential or about the scale factor behavior. These are obtained as outputs of the assumed linear relationship between the Hubble expansion parameter and the time derivative of the scalar field, once one integrates the field equations explicitely. As a consequence only a class of exponential potentials and their combinations can be treated. However, this is not a serious drawback of the method. In effect, combinations of exponentials are interesting alternatives since these arise in more fundamental (particle) contexts: supergravity and superstring, where these types of potentials appear after dimensional reduction. The relevance of the solutions found for the description of the cosmic evolution will be discussed in some detail making emphasis in the possibility for superquintessence solutions to arise. This issue has been alredy studied numerically in Ref..

We will be concerned with flat Friedmann-Robertson-Walker (FRW) cosmologies with the line element given by:

$$ds^2=dt^2+a(t)^2\delta _{ik}dx^idx^k,$$
(1)
where the indexes $`i,k=1,2,3`$ and $`a(t)`$ is the scale factor. We use the system of units in which $`8\pi G=c=1`$.

## 2 The Method

The Brans-Dicke field equations, in the presence of a self-interacting scalar field, are:

$$3H^2+3H\dot{\phi }\frac{\omega }{2}\dot{\varphi }^2=e^\phi (\rho _m+V),$$
(2)
$$2\dot{H}+3H^2+\ddot{\phi }+\dot{\phi }^2+2H\dot{\phi }+\frac{\omega }{2}\dot{\phi }^2=e^\phi \{(1\gamma )\}\rho _m+V,$$
(3)
$$\ddot{\phi }+\dot{\phi }^2+3H\dot{\phi }=\frac{e^\phi }{2\omega +3}\{(43\gamma )\rho _m+4V2V^{}\},$$
(4)
where we have introduced a new field variable $`\phi `$ instead of the original $`\varphi =e^\phi `$. In these equations $`\omega `$ is the Brans-Dicke coupling parameter, $`\gamma `$ is the barotropic index of the fluid of ordinary matter, $`V`$ is the self-interaction potential and $`H=\frac{\dot{a}}{a}`$ is the Hubble expansion parameter. The dot accounts for derivative in respect to the cosmic time $`t`$ meanwhile the comma denotes derivative in respect to the scalar field $`\phi `$. The energy density of the ordinary matter (cold dark matter plus baryons and/or radiation) is related with the scale factor through $`\rho _m=\rho _{0,\gamma }a^{3\gamma }`$, where $`\rho _{0,\gamma }`$ is an integration constant coming from integrating the conservation equation.

Let us introduce a new time variable $`dt=e^{\phi /2}d\tau `$. Then the field equations(2-4) can be rewritten in the following way:

$$3\overline{H}^2=\rho _m+\rho _\phi ,$$
(5)
$$2\dot{\overline{H}}+3\overline{H}^2=P_mP_\phi ,$$
(6)
$$\ddot{\phi }+\frac{1}{2}\dot{\phi }^2+3\overline{H}\dot{\phi }=\frac{43\gamma }{2\omega +3}\rho _m+\frac{2}{2\omega +3}(2VV^{}),$$
(7)
where

$`\rho _\phi ={\displaystyle \frac{\omega }{2}}\dot{\phi }^23\overline{H}\dot{\phi }+V,`$
$`P_\phi =\ddot{\phi }+{\displaystyle \frac{\omega +1}{2}}\dot{\phi }^2+\overline{H}\dot{\phi }V,`$
$`P_m=(\gamma 1)\rho _m.`$ (8)
Now the dot means derivative in respect to the new time variable $`\tau `$ and $`\overline{H}=\frac{da/d\tau }{a}`$. Combining of the above equations yields

$$\dot{\overline{H}}+\frac{1}{2}\overline{H}\dot{\phi }+3\overline{H}^2=\frac{(2\gamma )\omega +1}{2\omega +3}\rho _m+\frac{2\omega +1}{2\omega +3}V+\frac{1}{2\omega +3}V^{}.$$
(9)
An implicit symmetry of the LHS of equations (7) and (9) is made explicit by the change $`\overline{H}k\dot{\phi }`$, where $`k`$ is a constant parameter. Therefore, if we assume the linear relationship

$$\overline{H}=k\dot{\phi },a=e^{k\phi },$$
(10)
then, the following differential equation for the self-interaction potential $`V`$ is obtained:

$$V^{}+\frac{2\omega +14k}{2k+1}V=\frac{(43\gamma )k(2\gamma )\omega 1}{2k+1}\rho _m.$$
(11)
Straightforward integration of Eq. (11) yields

$$V(\phi (a))=Aa^n+Ba^{3\gamma },$$
(12)
where

$`n`$ $`=`$ $`{\displaystyle \frac{4k2\omega 1}{(2k+1)k}},`$
$`A={\displaystyle \frac{C_0}{3\gamma k\alpha }}`$ , $`B={\displaystyle \frac{\beta \rho _{0,\gamma }}{3\gamma k\alpha }},`$
$`\alpha =kn,\beta `$ $`=`$ $`{\displaystyle \frac{(43\gamma )k(2\gamma )\omega 1}{2k+1}}.`$ (13)
An interesting feature of this potential is that it depends on the type of ordinary fluid which fills the universe. Otherwise, it depends on the barotropic index $`\gamma `$ of the matter fluid. This fact implies some kind of interaction between the ordinary matter and the quintessence field much like the interacting quintessence studied in Ref..

By substituting Eq. (12) back into the Friedmann equation (5) one obtains

$$\overline{H}^2=\overline{A}a^n+\overline{B}a^{3\gamma },$$
(14)
where

$$\overline{A}=\frac{2k^2A}{6k^2+6k\omega },\overline{B}=\frac{2k^2(\rho _{0,\gamma }+B)}{6k^2+6k\omega }.$$
(15)
Finally we make another change of time variable: $`dr=a^{n/2}d\tau `$ to get Eq. (14) integrated in quadratures:

$$\frac{a^{\frac{3\gamma +n}{2}1}da}{\sqrt{\overline{B}+\overline{A}a^{3\gamma +n}}}=r+r_0,$$
(16)
where $`r_0`$ is another integration constant.

## 3 The Class of Solutions

Straightforward integration of Eq. (16) leads to:

$$a(r)=a_0\{\mathrm{sinh}[\mu (r+r_0)]\}^{\frac{2}{3\gamma +n}},$$
(17)
with

$$a_0=(\overline{B}/\overline{A})^{\frac{1}{3\gamma +n}},\mu =\frac{3\gamma +n}{2}\sqrt{\overline{A}}.$$
(18)
Then, from the linear relationship (10), one can derive the form of the scalar field:

$$\phi (r)=\frac{1}{k}\mathrm{ln}a_0+\frac{2}{(3\gamma +n)k}\mathrm{ln}\mathrm{sinh}[\mu (r+r_0)].$$
(19)
For purpose of observational testing we give the list of magnitudes of observational interest. The Hubble expansion parameter;

$$\overline{H}(r)=\frac{2\mu }{3\gamma +n}a(r)^{n/2}\mathrm{coth}[\mu (r+r_0)],$$
(20)
the matter density parameter,

$$\mathrm{\Omega }_m(r)=\frac{\rho _m}{3\overline{H}^2}=(\frac{3\gamma +n}{2\mu })^2\frac{\rho _{0,\gamma }}{3}\mathrm{cosh}^2[\mu (r+r_0)],$$
(21)
the scalar field density parameter,

$$\mathrm{\Omega }_\phi (r)=1\mathrm{\Omega }_m(r),$$
(22)
the deceleration parameter,

$$q(r)=1+\frac{1nk}{2k}+\frac{3\gamma +n}{2}\mathrm{\Omega }_m(r),$$
(23)
and the standard (perhaps not precisely the physically meaningful) equation of state of the scalar field,

$$\omega _\phi =\frac{P_\phi }{\rho _\phi }=1n/3.$$
(24)
This class of solutions depends (in principle) on six parameters: $`\gamma `$, $`\omega `$, $`k`$, $`\rho _{0,\gamma }`$, $`C_0`$ and $`r_0`$.

## 4 Reducing the Space of Parameters

A very crude analysis shows that $`\mathrm{\Omega }_\phi `$ should never be negative since, otherwise, this would imply coarse violations of the known energy conditions. Therefore $`\mathrm{\Omega }_m`$ should never be greater than unity. Since $`\mathrm{\Omega }_m`$ (Eq. (21)) is a maximum at $`r=r_0`$ where

$$\mathrm{\Omega }_m(r_0)=(\frac{3\gamma +n}{2\mu })^2\frac{\rho _{0,\gamma }}{3},$$
(25)
then by requiring that

$$\mathrm{\Omega }_m(r_0)=1,\frac{3\gamma +n}{2\mu }=\sqrt{\frac{3}{\rho _{0,\gamma }}}.$$
(26)
If we apply the normalization in which $`\rho _{0,\gamma }=1`$, then $`\frac{3\gamma +n}{2\mu }=\sqrt{3}`$ and Eq. (26) implies:

$$C_0=\frac{(6\gamma k^2+(4+3\gamma )k2\omega 1)(6k^2+6k\omega )}{6k^2(2k+1)}.$$
(27)
By considering Eq. (26) and evaluating Eq. (21) at present ($`r=0`$) one obtains:

$$r_0=\frac{1}{\mu }arccosh[\frac{1}{\sqrt{\mathrm{\Omega }_m(0)}}],$$
(28)
where, from observations, one may place $`\mathrm{\Omega }_m(0)=0.3`$.

After these considerations the space of parameters is reduced to the following free parameters: $`\gamma `$, $`\omega `$ and $`k`$. Besides the parameter $`\gamma `$ is fixed by the type of fluid dominating the given stage of the cosmic evolution. At present ($`r=0`$) $`\gamma =1`$, i. e., the universe is filled with a cold dust of dark matter plus baryons mainly.

We want to point out that all of the physical parameters of observational interest (equations (21-24)) can be given as functions of the redshift $`z`$ simply by replacing in the former equations

$$\mathrm{\Omega }_m(r)\mathrm{\Omega }_m(z)=\frac{(1+z)^{3\gamma +n}}{(1+z)^{3\gamma +n}+\mathrm{sinh}^2[\mu r_0]}.$$
(29)
This form of writing will be useful for purpose of observational testing.

## 5 Observational Testing

The main observational facts we consider are the following:

1.- At present ($`r=0`$) the expansion is accelerated ($`q(0)<0`$).

2.- The accelerated expansion is a relatively recent phenomenon. Observations point to a decelerated phase of the cosmic evolution at redshift $`z=1.7`$. There is agreement in that transition from decelerated into accelerated expansion occurred at a $`z0.5`$.

3.- The equation of state for the scalar field at present $`\omega _\phi 1`$ (it behaves like a cosmological constant).

4.- Although, at present, both the scalar (quintessence) field and the ordinary matter have similar contributions in the energy content of the universe ($`\mathrm{\Omega }_m(0)=0.3\mathrm{\Omega }_\phi (0)=0.7`$), in the past, the ordinary matter dominated the cosmic evolution,222The quintessence field $`\phi `$ should be subdominant in the past as required by nucleosynthesis constraints meanwhile, in the future, the quintessence field will dominate (it already dominates) and will, consequently, determine the destiny of the cosmic evolution.

In fig.1 the evolution of both the matter and scalar field density parameters is shown for $`\omega =1000`$ and an arbitrary $`k=530`$ ($`\gamma =1`$), as function of the redshift $`z`$. We see that the evolution of both magnitudes meets the observational requirements o matter dominance (scalar field subdominance) at high redshift (in the past) and scalar field dominance at low redshift (at present and in the future).

From fig.2, where the evolution of the deceleration parameter $`q\&z`$ is shown for the same $`\omega =1000`$ and $`k=530`$, we see that, at present ($`z=0`$), $`q(0)=0.55<0`$ (accelerated expansion) meanwhile, in the past ($`z>0`$) we had a stage of decelerated expansion. The transition from decelerated into accelerated expansion occurred at $`z=0.67`$, which is in acceptable agreement with the claimed value $`z0.5`$. Besides, the observational fact of decelerated expansion at redshift $`z=1.7`$ is fulfiled.

The equation of state is always a constant $`\omega _\phi =1.00007`$ for the former values of the free parameters $`\omega `$ and $`k`$, in good agreement with the accepted value. In general the solution found meets all of the above observational requirements (points 1-4) for a really wide range in the parameter space ($`\omega ,k`$).

## 6 Constrains from Superquintessence

As pointed out in Ref., superquintessence behavior ($`\omega _\phi <1`$) is a feature of models of quintessence based on non-minimally coupled scalar-tensor theories including Brans-Dicke theory with a self-interacting scalar field. In Ref., however, it was shown that expected superacceleration effect is small enough as to escape the reach of present and future observations. We now proceed to analyze the question of superquintessence within the frame of the solution found.

We shall interested in $`n0`$ (Eq. (13)) in order to have superquintessence behavior (see Eq. (24)). This, in turn, implies $`k>(2\omega +1)/4`$. Since, from experiment $`\omega 500`$, then $`k>250.25`$ is a necessary condition to have superquintessence. Consider $`n`$ as function of $`k`$ ($`\omega `$ constrained by experiment):

$$n=\frac{4k2\omega 1}{k(2k+1)},$$
(30)
which is a maximum $`n`$ at $`k=(2\omega +1+\sqrt{(2\omega +1)(2\omega +3)})/4`$, where

$$n=\frac{2\sqrt{(2\omega +1)(2\omega +3)}}{\{\frac{2\omega +1}{2}(1+\sqrt{\frac{2\omega +3}{2\omega +1}})\}^2+(\frac{2\omega +1}{2})(1+\sqrt{\frac{2\omega +3}{2\omega +1}})}.$$
(31)
For instance, for $`\omega =500`$, the maximum value for $`n`$ is $`n0.002`$ (at $`k500.75`$) so, according to Eq. (24), the maximum expected superacceleration effect is of very small magnitude ($`0.2\%`$). For other values of $`k>(2\omega +1)/4`$, this effect would be even smaller. However, as discussed in details in Ref., the physically meaningful definition of the equation of state for the scalar field reads;

$$\stackrel{~}{\omega }_\phi =1\frac{1}{3}\frac{\mathrm{ln}[\stackrel{~}{\rho }_\phi (r)/\stackrel{~}{\rho }_\phi (0)]}{\mathrm{ln}a(r)},$$
(32)
instead of Eq. (24). In (32) the following redefinition for the scalar field energy density has been used:

$$\stackrel{~}{\rho }_\phi =\rho _\phi 3\overline{H}^2\{1a_0^{1/k}(\mathrm{sinh}[\mu (r+r_0)])^{\frac{2}{(3\gamma +n)k}}\}.$$
(33)
In fig.3 we show the behavior of the equation of state both, the usual (Eq. (24)) and the physically meaningful (Eq. (32)), as function of the parameter $`k`$ for $`\omega =1000`$. Maximum expected superquintessence, in this case, is obtained for $`k800`$ (the standard equation of state) or $`k1300`$ (the physically meaningful equation of state). In both cases the effect is very small (of about 0.04$`\%`$).

## 7 Conclusions

We have found a new parametric class of exact cosmological solutions to Brans-Dicke theory of gravity with a self-interacting scalar field. A class of self-interaction potentials that is a combination of exponentials and that is dependent on the kind of matter fluid that fills the space-time arise. This potential suggests some kind of interaction between the quintessence field and the ordinary matter similar to the kind in Ref. .

The flat universe solution found meets the observational constrains on the physical parameters, among them $`\mathrm{\Omega }_m`$ ($`\mathrm{\Omega }_\phi `$), $`\omega _\phi `$, $`q`$ and the experimental requirement that $`\omega >500`$. It reflects the observational fact of a past decelerated, ordinary matter (including radiation) dominated universe and a future scalar field (dark energy) dominated universe.

The analytic study of the class of solutions derived in the present paper confirm the results of Ref. . Superquintessence behavior, although present for a given region in parameter space ($`\omega ,k`$), is behind the reach of present and future observational possibilities.

Acknowledgments

We thank our colleague Tame Gonzales for helpful comments and the MES of Cuba by financial support of this research.