Chapter 6
Solitons in Raman media

In this section, solitons supported by the Raman effect, rather than Kerr nonlinearity are considered. Raman scattering is the inelastic scattering of light, in which a photon is exchanged for a photon of slightly different energy. This energy exchange corresponds to a quantum transition within the medium, which is usually rotational or vibrational in nature. The Raman effect occurs when this transition is forbidden, but may happen indirectly via higher energy "virtual" levels. In the first stage of the process, a molecule accepts an incoming photon, and is excited to one of the virtual levels. In the second stage, the molecule decays from the virtual level to the second level, with the emission of the outgoing photon. This process can be subdivided into three separate phenomena (as shown in figure 6.1):

• Stokes Scattering occurs when the final state has a higher energy than the initial state. Therefore, the photons are effectively redshifted, and the medium absorbs energy.
• Anti-Stokes Scattering occurs when the final state has a lower energy than the initial state. Therefore, the photons are blueshifted, and the medium releases energy. Due to conservation of energy, anti-Stokes scattering must occur in tandem with another process that excites the medium. (Stokes scattering is one such process.) As the molecules can revert to the ground state by other means, anti-Stokes scattering will usually occur to a lesser extent than Stokes scattering, and so the ultimate effect of Raman scattering is usually a net redshifting.
• Rayleigh Scattering occurs when the final state has the same energy as the initial state. Therefore, no energy is transferred, and hence Raman scattering has not occurred. This in fact is the most common outcome, but it results in no change to the photon energy, and so will not be considered further.

The Raman effect requires the quantum coherence between the two energy levels to have an appropriate value. (This quantity will be defined in section 6.1, when the equations of motion are derived.) In this chapter, Stimulated Raman Scattering (SRS) is considered, whereby the coherence is provided by light already at the Stokes or anti-Stokes frequency. The Raman effect will then cause this frequency to be amplified, at the expense of the pump frequency. (Not surprisingly, optical amplifiers based on this principle are widespread [8].) In addition to SRS, the Raman effect can be initiated by vacuum noise, in what is known as Spontaneous Raman Scattering. A Stokes line or anti-Stokes line can itself undergo Raman scattering, to create a further line. These lines can create more lines, and so on, generating a comb of frequencies in what is known as cascading Raman scattering.

Like Kerr nonlinearity, SRS can support a variety of soliton solutions. These include hybrid bright-dark temporal solitons in which a bright-soliton of Stokes-shifted light co-moves with a dark-soliton in the optical pump [38], and more conventional bright solitons [39, 216, 217]. It is also possible to see solitons where the Raman effect and the Kerr effect are combined [218], but here we consider pure Raman solitons.

A Raman soliton consisting of three or more discrete frequency components requires that these components be phase matched. When phase mismatch is present, instability will occur, and the soliton will break up [219]. In this chapter, however, a new class of solutions are found, which can exist in the presence of strong phase mismatch.

The solutions derived aren’t applicable to silicon waveguides. This is due to the fact that they require the absence of Kerr nonlinearity, which is not the case for silicon. Furthermore, the Raman effect in silicon is difficult to observe; it has a very narrow bandwidth, and only becomes important in devices designed to exploit it [76, 77, 78, 79], where a signal frequency precisely tuned to the Raman-shifted pump frequency is amplified. In fact, from the point of view of solitons, the Raman effect in silicon is largely irrelevant [116], which is why it was not included in the model derived in chapter 2. When it is considered (e.g. [60]), it is usually represented using a modified Kerr nonlinearity term with hysteresis [60], rather that the coupled envelope model derived below. Finally, the solitons require an exotic state consisting of two levels with equal populations in each, which would not be possible in silicon.

A more suitable medium is hydrogen or deuterium [220, 221]. In addition to being placed in a gas cell, these gases can be placed in a hollow-core photonic crystal fibre [222]. These fibres have a hole running along their entire length, surrounded by a structured region of glass and air spaces, designed so that the majority of the light travels through the hole. Using such an arrangement, the gas can provide the Raman activity, whilst the fibre provides a tightly controlled dispersion relation.

6.1 Equations of stimulated Raman scattering

We now need to derive a set of equations describing how light evolves in a Raman active medium. (This derivation is shown in [97], but only in the much simpler case where only two field components are present.) We will assume that there are two real energy levels in the medium, plus N – 1 Raman levels which mediate scattering between N frequency components in the fields. This is shown in figure 6.2.

For clarity, we will use upper-case omega (Ω) to represent the frequencies corresponding to quantum energy levels, and lower case omega (ω) to represent the frequencies of components in the electric field. We number these frequencies starting from the highest (i.e. Ω₁, ω₁) and decreasing to the lowest (i.e. Ω_N–1, ω_N). At no point do we make a distinction between the numbering of Stokes, anti-Stokes, or pump frequencies. Therefore we can arbitrarily choose one of frequency components to be the pump frequency, in which case all of the components with lower ordinal numbers become anti-Stokes lines, and those with higher ordinal numbers become Stokes lines.

As we know the allowed states of the system, it is convenient to use a matrix formulation of quantum mechanics. We start by assuming that the states of the system, plus the virtual states form a complete orthonormal basis set. Therefore we can write an arbitrary wavefunction |ψ⟩ in this basis as

(6.1)

where |g⟩ is the ground state, |e⟩ is the excited state, and |n⟩ (where n = 1,2,,N – 1) is the nth Raman level. The coefficients c_g, c_e and c_n correspond to the relative amplitude of each state. This mixed state can be represented by a column vector

(6.2)

As we are using the eigenvectors of the unperturbed Hamiltonian as our basis, the matrix form of the unperturbed Hamiltonian will be a diagonal matrix consisting of the corresponding eigenvalues (i.e. the energies of those states). Therefore, we obtain

(6.3)

where we have defined the ground state to have zero energy, ℏΩ_R is the energy of the excited state, and ℏΩ_n is the energy of the nth Raman level. The frequency Ω_R is the frequency change caused by the Raman scattering.

The interaction Hamiltonian for a dipole transition from state x to state y is –⟨y| ⋅|x⟩ where is the dipole moment associated with that transition. Therefore, if we assume that each allowed transition is a dipole transition with a dipole moment of μ (whilst all the forbidden transitions have a dipole moment of zero), then the interaction Hamiltonian is

(6.4)

where E is the electric field, and is the polarisation operator, which is given by

(6.5)

Assuming that the number of Raman-active sites per unit volume (N₀) is very large, then the nonlinear part of the polarisation density (as defined in section 2.2.3) can be given statistically as

(6.6)

where ⟨⟩ is the expectation value of the polarisation operator. This can be calculated in terms of the density matrix , giving

(6.7)

Substituting the individual elements of and into this gives

(6.8)

where subscript g and e correspond to the ground and excited states respectively, and subscript n corresponds to the nth Raman level. The evolution of a density matrix is given by

(6.9)

where Ĥ is the full Hamiltonian (given by Ĥ₀ + Ĥ_int). We now need to reduce this into equations describing how the electric field evolves.

Evolution of the electric field

We proceed by taking two of the (generalised) components of equation 6.9, specifically

(6.12)

where A_n is the slowly varying envelope for the nth line, with frequency ω_0n and wavenumber β_0n. Instead of directly substituting this into equations 6.10 and 6.11, we will first prepare for what is known as the Rotating Wave Approximation. We start by representing the density matrix elements in the form

(6.13)

where ρ'_xy is the density matrix in a Rotating Wave Basis, Ω_x and Ω_y are the De Broglie frequencies of the given levels, and K_x and K_y are the De Broglie wavenumbers. (In the above notation, Ω_g = 0 and Ω_e = Ω_R). The reason for this is as follows: If there is no interaction between two levels, then their coherence will be precisely e^{i(K_x–K_y)z–i(Ω_x–Ω_y)t}, which is the beat oscillation between the two wavefunctions. If the levels do interact, then this oscillation will remain, but it will be modulated by ρ'_xy, which in general is slowly varying. Therefore, we can split the coherence into a rapidly oscillating part (e^{i(K_x–K_y)z–i(Ω_x–Ω_y)t}) which accounts for the innate quantum oscillation between the wavefunctions, and a slowly varying part ρ'_xy which accounts for non-trivial interactions between the levels.

Transforming into this basis, and collecting all the oscillating terms onto one side yields

(6.24)

gives the coherence between the energy levels, as was mentioned above. The ρ_ee and ρ_gg terms in equations 6.22 and 6.23 correspond to an oscillating Stark shifting of the levels [97]. (The Stark effect is the phenomenon whereby the energy levels of a quantum system are shifted by an applied electric field.) This term can be removed, as it oscillates with the electric field, and thus will integrate to zero over a short time span [97]. We can then substitute the remaining terms of the matrix elements into equation 6.8 to give

(6.25)

Assuming that the dipoles point in the same direction as the electric field (specified by '), this can be rewritten as

(6.26)

where μ is the dipole magnitude. This is very similar in structure to the polarisation in slowly varying envelope form, which (from equation 2.60 in section 2.2.3) is

(6.27)

Therefore, by matching together the parts which oscillate at the same frequency, we can identify the polarisation envelopes as

(6.28)

Substituting these envelopes into equation 2.58 gives

(6.29)

where β_mn is the mth dispersion coefficient for the nth field envelope. It should be noted that the summation index now starts from 0 rather than 2, as the β_0n and β_1n coefficients will in general be different.

We have also replaced the ω_0n coefficients with a single value of ω₀, to remove the explicit frequency dependence of the nonlinearity. This is valid so long as the spectral range of the frequency comb is relatively narrow. (In principal we could have a formulation where the nonlinearity terms were different for each line. However, to do this properly, it would also be necessary to replace the dipole magnitude μ and the detuning δ with a separate value for each possible transition. This would require a substantial amount of physical data, and would also greatly complicate mathematical analysis. It is sufficient for our purposes to replace the nonlinearity coefficient with a single value that can be scaled away.)

As the individual spectral peaks span a frequency range that is narrower than the whole comb, the effect of group velocity dispersion (GVD) will predominantly occur between the spectral peaks, rather than within them. We can therefore make a common [218], but not universal [217], approximation of removing the terms in β_2n to give

(6.30)

This does not amount to ignoring the GVD, as it still manifests itself in the differences between the values of β_1n.

Evolution of the Raman coherence

It is now necessary to determine how the Raman Coherence evolves. This can be done by considering the ρ_eg component of equation 6.9, which is

(6.31)

Replacing ρ_eg with Q (using equation 6.24) gives

(6.32)

Substituting in the values for ρ_ng and ρ_ne (using equations 6.22 and 6.23), and resolving E into a set of envelope functions (using equation 6.12) gives

(6.34)

The Q|A_n+1|² – Q|A_n|² terms sum to zero, and hence we are left with

(6.35)

where σ (defined as σ ≡ ρ_ee – ρ_gg) is the Population Inversion Parameter.

Evolution of the population inversion parameter

It is now necessary to determine how the population inversion parameter evolves. This can be done by considering the ρ_ee and ρ_gg components of equation 6.9, which subtracted from one another give

(6.36)

By an essentially identical process to that used in deriving the equation for Q (i.e substituting in the equations for the electric field and matrix elements, and then selecting the resonant terms), we obtain

(6.37)

This assumes that deexcitation can only happen via anti-Stokes scattering. In general, however, other processes can deexcite the medium, and so a term of the form r should be added to the right hand side of equation 6.37, where r is the decay rate.

Conservation of population

It should be noted that the quantity σ² + |Q|² is an integral of motion. If we decompose it into the constituent wavefunctions, we obtain

(6.38)

where |g⟩ and |e⟩ are the wavefunctions of the ground and excited states of the system. Therefore, the integral of motion corresponds to the conservation of the population of states. The wavefunctions are normalised, and so we require that

(6.39)

Equations in intrinsic units

The equations for N-field Raman scattering (equations 6.30, 6.35 and 6.37) can be made dimensionless by choosing intrinsic units such that

An example of a pure Raman medium (i.e. one that conforms to the above model, without complications due to Kerr nonlinearity and other physical effects) is deuterium gas. At a pump frequency of ω_n = 3.71 × 10¹⁴, this has a Raman coupling coefficient of R = 4.5 × 10^–7 C²m²J^–2s^–1 [220]. At atmospheric pressure (where the volumetric density N₀ = 2.45 × 10²⁵ m^–3 m^–3), this gives a distance unit of Z₀ = 1mm.

6.2 Raman solitons with phase mismatch

When a Raman soliton consists of only two spectral lines, the phase mismatch is irrelevant, as it can be scaled away by switching to a moving frame of reference. With three or more lines, however, it is impossible to remove phase mismatch in the general case. Existing multi-line solutions simply ignore the phase mismatch parameters [223, 224], and so are only valid in the limit of short pulses. (The reason for this limit can be seen by noting that the phase mismatch parameter is in fact a measure of phase mismatch per unit length. Therefore, the total phase shift in a pulse is given by θ_n = η_nΔζ, where Δζ is its approximate spatial length. For the phase shift to be irrelevant, it must be much less than a wavecycle, giving the condition θ_n ≪ 2π, and thus Δζ ≪ 2π/η_n.)

In this section, these solutions are extended to include dispersion relations with phase mismatch.

6.2.1 Band-gaps and tail analysis

Before attempting to find soliton solutions, we will explore the preconditions for soliton existence. Firstly, we need to know the linear dispersion relation (and the position of the band-gap that exists within it), so that we can avoid the instability resulting from an intersection between a soliton’s dispersion relation and a linear mode.

The dispersion relation of linear waves can be obtained by substituting B_n = b_ne^ikζ–iωτ into equations 6.45. The complex constants b_n give the amplitudes and phases of each field component. This gives

(6.48)

where the matrix M_N is defined as

(6.49)

Therefore, the dispersion relation is given by |M_N| = 0, which can be written recursively as

(6.50)

where |M₁| = η₁ – k + v₁ω and |M₀| = 1. This equation will trace out N curved paths in space, which correspond to the possible linear wave dispersion relations. For very large values of ω or k (in both the negative and positive directions), the equation approaches the limit

(6.51)

and so the set of velocities (ω/k) for linear waves in this limit is

(6.52)

Therefore, at each extreme of k or ω, the linear wave dispersion curve velocities must approach a set of values that matches up to the set V . We can sort the linear wave dispersion curves given by equation 6.50 into three distinct categories, by using this knowledge of how the curves behave at their extremal values. Each individual curve can:

1. Change from a positive velocity at ω = –∞ to a negative velocity at ω = ∞, reaching a maximum value of k somewhere in between. (A negative velocity simply means a retardation with respect to the moving frame of reference, rather than backwards propagation.)
2. Change from a negative velocity at ω = –∞ to a positive velocity at ω = ∞, reaching a minimum value of k somewhere in between.
3. Have the same sign of velocity at both extremes of ω and so cover the entire range of k.

(This excludes the unimportant degenerate cases where v_n = v for one of the field components.)

Figure 6.3: Linear dispersion curves for a variety of N. For even N, a bandgap exists, whilst for odd N, the bandgap is obstructed. Shown for ηn = 0 and vn = n – –.

When all of the curves fall into the first and second categories, there will exist a band-gap between the maxima and minima. This allows a soliton to exist without its dispersion relation overlapping with that of a linear wave mode. When a bandgap does exist, it allows for a range of soliton velocities between the least negative and the least positive values of v_n.

It should be noted that for this band gap to exist, exactly half of the v_n coefficients must be positive, and the other half negative, so that each of the dispersion curves can have a negative gradient at one end, and a positive gradient at the other. (As we are working in a moving frame of reference, a negative group velocity doesn’t necessarily correspond to a negative group velocity in the lab frame.) This precludes the existence of a band gap when the number of field components is odd, as at least one dispersion curve must fit into the third category. In figure 6.3, linear wave dispersion relations for a selection of N field systems are shown, where η_n = 0 and v_n = n – N/2 – 1/4. For even N a bandgap exists, whilst for odd N it is obstructed.

Tail analysis

If we assume that soliton solutions have exponentially decaying tails, we can investigate these tails by substituting

(6.53)

into the equations of motion, where τ'≡ τ – v_sζ is the retarded time for a soliton moving at (reciprocal) speed v_s, and λ is the exponential constant, which should have a positive real part for the leading edge of the pulse, and a negative real part for the trailing edge. It should be noted that this procedure is very similar to the above bandgap analysis, because mathematically speaking, an exponential decay is equivalent to a sinusoidal oscillation with imaginary frequency. Therefore, we can derive the same equation by the making the transformation

(6.57)

where the matrix L_N is defined as

(6.58)

The solutions for λ are given by |L_N| = 0, which can be written recursively as

(6.59)

where |L₁| = η₁ –κ–iλ and |L₀| = 1. This is an Nth degree polynomial in λ, and so λ has N solutions. The polynomial is unchanged by the transformation λ →–λ^∗. Therefore, solutions must exist in pairs of the form . When N is even, each value of λ is generally complex. However, when N is odd, there must be a lone solution of λ that pairs up with itself. To satisfy the relationship between pairs, this root must lie on the imaginary-axis. A purely imaginary decay constant corresponds to a simple harmonic wave (or in the degenerate case of λ = 0, a constant amplitude). As with the band-gap analysis, this suggests instability for odd N.

6.2.2 Analytical soliton solutions

We can obtain analytical solutions by making the following ansatzes:

• Q can be written in the form Q = α + iq, where q is a real function, and α is a real constant.
• The fields can be written in the form

  (6.60)

  where f is a function common to all the fields, x ≡ c is the retarded time scaled by a width parameter c, and κ is the wavenumber. The amplitudes and phases of each soliton are given by the complex constants b₁,b₂,b_N.

Substituting these into equations 6.45, 6.46 and 6.47 gives

(6.64)

By making the further ansatz that

(6.65)

we can manipulate equation 6.61 into the form

(6.66)

where the matrices S_N and T_N are given by

As S_N and T_N are both real, all the values of b_n will have the same complex phase, and so γ will be real. Therefore, from equations 6.62, 6.63 and 6.65 we obtain the simplified equations

(6.74)

In this form, the population of the system is given by

(6.75)

and hence the population normalisation condition is given by

(6.76)

Substituting equations 6.74 and 6.76 into equation 6.73 yields an uncoupled equation for σ

(6.77)

This can be separated to give

(6.78)

Evaluating these integrals yields

(6.79)

where C is the integration constant. Rearranging this to give an equation in σ is not possible in closed form. However, if we assume that J = 0, the equation reduces to

(6.80)

Physically speaking, the assumption of J = 0 means a pulse of light is being fired into a medium with an equal population between the two states, and perfect coherence between those states. It is possible to prepare such a medium [224], and so this is a realistic scenario. The equation can now be solved to give

(6.81)

where we have set C = 0, as this merely shifts the solution in time. The solution can be written as

(6.82)

where

(6.83)

This in turn can be written as

(6.84)

where the F_n function is defined as

(6.85)

This function is a generalisation of the familiar sech function, to which it is identical when n = 1. Substituting equation 6.84 into equation 6.72 gives an uncoupled equation for q

(6.86)

This has a solution

(6.87)

and so the coherence can be written as

(6.88)

Substituting equation 6.84 into equation 6.74 gives a solution for |f|². The equations of motion are satisfied by taking the real positive route of this function, to give

(6.89)

Before solving equations 6.69 and 6.70 (in order to calculate the actual field envelopes) we should note that it is homogenous, and so the values b_n will have an arbitrary multiplicative constant, m. From its definition, γ is proportional to m², and so f is proportional to m^–1. Therefore, when we multiply f by b_n to obtain the field envelopes, m will cancel out. We can remove m from our calculations entirely by specifying that the values of b_n must satisfy the condition

(6.90)

We can then write the field envelopes as

(6.91)

where the coefficients b_n have yet to be found.

This type of solution is notable for carrying a topological charge [223]. The sech-like Kerr soliton considered in section 3.1 is non-topological, as it can be continuously deformed into a non-solitonic (uniformly zero) solution, by letting the wavenumber tend towards zero. The Raman soliton, however, cannot be deformed like this, as the coherence Q is different between the extremes of time, and no form of scaling will remove it.

6.2.3 Calculation of field amplitudes

Equations 6.69 and 6.70 impose contradictory requirements on the values of b_n. However, a limited set of solutions can be obtained by solving equation 6.69 (as is shown below), and then matching T_N to this solution by explicitly requiring that the phase-mismatch coefficients conform to

(6.92)

The constants α and κ can be chosen arbitrarily (subject to the above condition that –1 < α < 1), and hence can be used to gain a certain amount of control over the phase mismatch parameters. It should be noted that in the case α = κ = 0, the phase mismatch coefficients all become zero. This special case is already known in the literature [224], but the phase mismatched generalisation is novel.

The reciprocal group velocity parameters were chosen as

(6.93)

where the parameter ξ is the group velocity dispersion, with positive values giving normal dispersion, and negative values giving anomalous dispersion. The values of v_n are shifted in order to guide the soliton’s spectrum through the bandgap derived in section 6.2.1. (In the case of odd N, when the band gap is always occluded, only one linear band will intersect with the soliton spectrum.) The constant ϕ can lie within the range – < ϕ < .

If 6.69 is to have non-zero solutions in b_n, then the soliton width c, must be chosen such that.

(6.94)

It is helpful to consider the soliton velocity in relation to the group velocity dispersion, and so we define the proportional soliton velocity v_p, by the identity v_s ≡ ξv_p. This enables us to write S_N as

(6.95)

Therefore, we can determine the soliton width for all values of ξ by solving

(6.96)

for ξc. It is apparent that the width, 1/c, of a soliton is proportional to the group velocity dispersion parameter. Furthermore, the values of b_n are independent of ξ. Solutions for ξc were found numerically. For certain parameter sets, the solutions are all non-real, and so the soliton solutions are not applicable. The regions where real solutions for c exist are shown in figure 6.4.

For cases where real values of c exist, the values of b_n (subject to condition 6.90) were calculated numerically. Solutions at v_s = 0.05, for the six first applicable cases are shown

Substituting these values of b_n and c into 6.91, gives a soliton solution. This solution is a generalisation of the solutions given in [224].

Gaussian profile of the field amplitude coefficients

As N increases, the shape of the b_n distribution tends towards a Gaussian profile. This can be shown by considering an arbitrary row of equation 6.69, such that

(6.97)

This can be rearranged to give

(6.98)

For a large number of fields, b_n can be approximated by a continuous function b(n), where the values of b are approximately equal for adjacent values of n. Similarly, v_n can be represented by a continuous function v(n). Therefore, we can approximate the left hand side of equation 6.98 with a derivative, giving

(6.99)

This can be solved to give

(6.100)

(where the term has been absorbed into the arbitrary constant of the remaining integral). For the values of v given by equation 6.93, we obtain

(6.101)

where m is the integration constant. This is a Gaussian function centred about the zero group velocity point (at n = N – ϕ) of the dispersion relation. Notably, to see this Gaussian distribution (and hence to have spectral localisation), the GVD parameter ξ must be positive. Therefore, unlike Kerr solitions, this class of soliton requires normal GVD.

Phase mismatch parameters

Once the values of b_n have been calculated, the permitted values of η_n can be trivially calculated from equation 6.92. It follows from the Gaussian profile of b_n that the phase mismatch parameters will have a profile of the form

(6.102)

where the approximation N ≫ 1 has been made. This is a hyperbolic cosine function centred about the zero-velocity point (at n = N/2 – ϕ) of the dispersion relation.

6.2.4 Simulation of soliton propagation

The full equations were solved numerically, using the soliton solutions as an initial condition at ζ = 0, and then integrating equations 6.45, 6.46 and 6.47 using the methods described in appendix A.4. Random noise with a relative amplitude of 1% was added to the initial fields, in order to demonstrate soliton stability under perturbation. The simulation results are displayed in terms of the magnitude of the electric field, assuming that the Raman frequency is Ω_R = 10. In each simulation, the parameters ϕ = –0.25 and v = 0.05 were chosen. As a proof of concept, the propagation of these solutions was tested with (essentially arbitrary) non-zero values of α and κ, and hence of η_n. A ten field soliton was chosen, thus avoiding the potential instabilities that can occur for an odd number.

With zero-phase mismatch (α = 0, κ = 0 in figure 6.8), the soliton propagates stably, despite the addition of noise. When phase match is introduced, slight perturbations appear in the time domain. Despite this (and despite the very long propagation distance) the pulse remains localised. Furthermore, in the frequency domain the amplitude of each field component is constant, indicating that no net Stokes or anti-Stokes shifting is occurring. Therefore, the waveforms are almost certainly solitons.