## Abstract

We consider the step Riemann problem for the system of equations describing the propagation of a coherent light beam in nematic liquid crystals, which is a general system describing nonlinear wave propagation in a number of different physical applications. While the equation governing the light beam is of defocusing nonlinear Schrödinger (NLS) equation type, the dispersive shock wave (DSW) generated from this initial condition has major differences from the standard DSW solution of the defocusing NLS equation. In particular, it is found that the DSW has positive polarity and generates resonant radiation which propagates ahead of it. Remarkably, the velocity of the lead soliton of the DSW is determined by the classical shock velocity. The solution for the radiative wavetrain is obtained using the Wentzel–Kramers–Brillouin approximation. It is shown that for sufficiently small initial jumps the nematic DSW is asymptotically governed by a Korteweg–de Vries equation with the fifth-order dispersion, which explicitly shows the resonance generating the radiation ahead of the DSW. The constructed asymptotic theory is shown to be in good agreement with the results of direct numerical simulations.

## 1. Introduction

Dispersive shock waves (DSWs), also termed undular bores in fluid mechanics, are generic solutions of nonlinear dispersive wave equations, including the Korteweg–de Vries (KdV), nonlinear Schrödinger (NLS) and Sine–Gordon equations. A DSW forms due to the dispersive resolution of a discontinuity and is the dispersive equivalent of a gas dynamic shock for which a discontinuity is resolved by viscosity [1]. A DSW is a non-steady modulated wavetrain which continually expands and has solitary waves at its leading edge and linear, small amplitude waves at its trailing edge (for the case of negative dispersion; if dispersion is positive then the orientation of the DSW, i.e. the relative positions of the linear and soliton edges, swap). This modulated wavetrain provides an oscillatory transition between the two levels of the initial discontinuity.

DSWs/undular bores are a common wave form which can be found in a broad array of physical systems. The classical undular bore is the tidal bore found in regions of large tidal flows and suitable topography, for example, the Severn Estuary in England and the Bay of Fundy in Canada. However, undular bores arise in a wide range of fluid systems, including the atmosphere, an example being morning glory clouds [2,3] and the semi-diurnal internal tide [4]. They also arise in geophysics (magma flow) [5–7] and Fermi gases [8]. Of particular relevance to the present work, they arise in nonlinear optics for a wide range of optical materials, including photorefractive crystals [9–11], optical fibres [12–17], nonlinear thermal optical media [18,19], colloidal media [20,21] and nematic liquid crystals [21,22].

DSW solutions of nonlinear dispersive wave equations are usually found using Whitham modulation theory [1,23,24]. Whitham modulation theory is a method for analysing slowly varying (modulated) wavetrains and deriving equations for the parameters, mean height, wavenumber, amplitude, etc., of such wavetrains. It is equivalent to the method of multiple scales, but much simpler than this to implement. When the underlying wavetrain is stable, the modulation equations form a hyperbolic system for the wavetrain parameters. It was found that a simple wave solution of the hyperbolic modulation equations for the KdV equation corresponds to a DSW [25,26], so that the standard method for finding DSW solutions is from the modulation equations for the relevant governing equation. This original method due to Gurevich & Pitaevskii [25] and Fornberg & Whitham [26] relies on the hyperbolic modulation equations being in Riemann invariant form, which is guaranteed if the governing equation is integrable with an inverse scattering solution [27]. However, most equations governing DSWs in physical applications are not integrable. This limitation was overcome to a certain extent when it was found that the leading, soliton, edge and trailing, small amplitude wave, edge of a DSW could be determined without a knowledge of the full Whitham modulation equations [28].

In this work, a DSW due to coherent light propagation in a nematic liquid crystal is analysed. While the specific context is light propagation in a nematic liquid crystal, equations similar to those for light propagation in this medium also arise for other nonlinear optical media [18,29–34], in fluid mechanics [35] and in models of quantum gravity [36]. An optical DSW in a nematic liquid crystal is found to possess a number of unique features. While the equation governing the optical field in a nematic liquid crystal is of defocusing NLS-type [37], the DSW is found to be of positive polarity, KdV-type, due to the effect of the nematic optical medium, which has a highly ‘non-local’ response [38–40]. It is further found that the dispersion relation for linear waves is non-convex, so that there is a resonance between the DSW and dispersive radiation. This results in a resonant wavetrain propagating ahead of the DSW. A similar resonant coupling between a DSW and radiation was found for nonlinear optical beam propagation in optical fibres when higher order dispersive terms were included in the governing NLS equation to enable such coupling, both without [14,15,17] and with [16] loss. The driving mechanism is the resonant coupling with higher order dispersion, which can also occur with just a soliton [41]. The total structure of the Riemann problem solution is then found to consist of four distinct regions, (i) an expansion wave linking the initial level behind to an intermediate shelf, (ii) a KdV-type DSW on this shelf, (iii) a resonant wavetrain leading the DSW, and (iv) a front bringing the resonant wavetrain down to the initial level ahead. Asymptotic solutions for all these four regions are obtained and compared with full numerical solutions of the governing equations, with generally excellent agreement being found.

The paper is organized as follows. In §2, the equations governing light beam propagation in a nematic liquid crystal are introduced and related to similar systems of equations in other physical contexts. In §3, the dispersive-hydrodynamic properties of these nematic equations are analysed and it is found that, while the dispersionless limit is described by a hyperbolic system equivalent to the shallow water equations, which is consistent with the dispersionless limit of the defocusing NLS equation, the linear dispersion relation is non-convex, implying the possibility of the formation of a KdV-type DSW in the low-frequency region and the generation of high-frequency resonant radiation by the DSW. This effect is a counterpart of the well-known radiating solitary waves in systems with higher order dispersion studied previously in many physical contexts, from gravity–capillary waves [42] to optical supercontinuum generation (see [43] and references therein). In §4, the fifth-order KdV equation (also known as the Kawahara equation) is derived from the nematic equations under a balance between strong non-locality and the small amplitude, long-wave approximation. The coefficient of the fifth-order dispersion term is proportional to the non-locality squared. It is then shown numerically that the effect of the non-locality on the DSW is the generation of a radiative wavetrain ahead of the DSW. In contrast to the well-studied radiating solitons of the fifth-order KdV equation, which are intrinsically unsteady, the solitary wave at the leading edge of the radiating DSW remains steady due to energy influx from the rest of the DSW. It can then be well approximated by the standard KdV soliton if the higher order dispersive term is sufficiently small. In contrast to previous work [22], it is found that the velocity of the leading edge of the KdV DSW is given by a classical shock jump condition, rather than the conservation of Riemann invariants [28]. This suggests that the resonant wavetrain acts as an effective viscous loss term for the DSW. In §5c, a Wentzel–Kramers–Brillouin (WKB) solution is constructed for the rapidly oscillating, resonant, linear radiative wavetrain in the full nematic system under the assumption that the lead solitary wave in the DSW can be approximated by a KdV soliton. Section 6 is devoted to comparisons of the constructed modulation solution with full numerical solutions of the nematic system.

## 2. Nematic equations

In this paper, we consider the propagation of a polarized, coherent beam of light through the medium of a nematic liquid crystal [39,40]. We assume that the electric field of the light is in the *x*-direction and that the beam propagates in the *z*-direction. Nematic molecules are elongated molecules, hence their name as nematic comes from the Greek word for thread, along which electrons can move freely. Hence an electric field, either an external static electric field or the electric field of light, results in the nematic molecules becoming dipoles and rotating in the direction of the electric field due to the resulting torque in order to minimize the potential energy [39,40]. The molecules rotate until the elastic forces balance the electrostatic forces. This rotation changes the refractive index of the nematic medium. Normally, a nematic is a focusing medium, so that the refractive index increases on rotation of the molecules. This self-focusing can then balance the diffractive spreading of a light beam, so that a bright optical solitary wave, termed a nematicon, can form [39,40,44]. However, the addition of azo-dyes to the nematic medium changes its structure so that it can become a defocusing medium as rotation of the molecules then decreases the refractive index [45]. In this case, a dark solitary wave, a dark nematicon, can form, a dip in a uniform background, rather than the rise from a background of a bright nematicon in the focusing case. The added complication of the nematic medium is that if the nematic molecules are initially aligned with their axis, termed the director, orthogonal to the electric field, the optical Freédericksz threshold exists so that a minimum electric field strength is required to overcome the elastic forces of the nematic medium before the molecules can rotate [46,39]. To enable nematicons to form at milliwatt power levels so that there is not excessive heating of the nematic, which can result in it undergoing a phase transition, an external static electric field is applied to pre-tilt the nematic molecules at an angle *θ*_{0} to the *z*-direction. In the particular case *θ*_{0}=*π*/4, the Freédericksz threshold vanishes [44,39].

Let us denote the extra rotation from the pre-tilt caused by the electric field of the light beam to be *θ*. Then in the paraxial, slowly varying envelope approximation the system of equations governing the propagation of a nonlinear light beam through a defocusing nematic liquid crystal is [38–40,45]
*u* is the complex valued envelope of the electric field of the light beam. The parameter *ν*, termed the non-locality, measures the elastic response of the nematic and is large, *ν*=*O*(100), in experiments [47]. This large value of the non-locality *ν* will be found to have a dominant effect on the structure of a DSW in a defocusing nematic liquid crystal. The parameter *q* is proportional to the square of the pre-tilting electric field. The electric field equation (2.1) is an NLS-type equation, which is coupled to equation (2.2) for the response of the nematic medium.

The context of the system of equations (2.1) and (2.2) has been explained in detail in terms of the nonlinear optics of liquid crystals. However, this system arises in a wide range of applications. In nonlinear optics, it arises whenever the response of the optical medium is based on some type of diffusive phenomenon [31], for example, it arises in the optics of nonlinear thermal media [30,18], for example, lead glasses [29,33,34] and certain photorefractive crystals [32]. A similar system of equations arises in simplified models of fluid turbulence [35] and quantum gravity [36].

In this paper, we consider the Riemann problem for the nematic system (2.1) and (2.2). The electric field equation (2.1) will be solved with the initial condition
*z*=0, with *u*_{3}>*u*_{1} so that a DSW is be generated. For consistency, the director equation (2.2) gives
*z*=0.

A typical solution of the nematic equations for the step initial condition (2.3) for large non-locality *ν* is displayed in figure 1*a*. For comparison, the solution for small *ν* is displayed in figure 1*b*, noting that the nematic equations (2.1) and (2.2) reduce to the NLS equation in the limit *ν*→0. As found in previous work [22], for large values of the non-locality *ν* the solution does not display the typical defocusing NLS DSW structure of figure 1*b* [37], even though the electric field equation (2.1) is of defocusing NLS-type. There is a KdV-type DSW in the electric field on the intermediate shelf of height *u*_{2} between the initial levels *u*_{3} and *u*_{1}. Preceding this DSW, there is a relatively high-frequency wavetrain, with a front which brings it back to the initial level *u*_{1}. The KdV DSW and resonant wavetrain are mirrored in the director response, at a much reduced amplitude, with the resonant wavetrain in the director having amplitude *O*(*ν*^{−1}). The inset in figure 1*a* shows the details of this resonant wavetrain. In this paper, the complex wave structure seen in figure 1*a* is understood as a radiating DSW. Such radiating DSWs typically arise for nonlinear wave equations with higher order dispersion, the model equation being the fifth-order KdV equation, or Kawahara equation. Although the theory of radiating solitons for the fifth-order KdV and similar equations is well understood (see [48–50] and references therein) the counterpart for DSW theory has only started to be explored (see the monograph [51] and references therein and the recent papers [12,14–17]). In addition to the leading resonant wavetrain, there are also radiative waves on the intermediate shelf on which the DSW sits. These are most likely due to internal resonances within the DSW which is a modulated periodic wave with a range of phase and group speeds.

## 3. Nematicon dispersive hydrodynamics

To analyse the Riemann problem (2.1)–(2.4), it is instructive to introduce the Madelung transformation
*ν*: the long scale *O*(*ν*^{1/2}) and the short scale *O*(1), which is consistent with the two distinct types of oscillatory structures observed in figure 1*a*. These distinct structures are characterized by differing typical wavelengths and different types of dispersion, which can be understood by analysing the linear dispersion relation for this system.

Linearizing the hydrodynamic form of the nematic equations (3.2)–(3.4) around the background levels *Θ*_{1} due to the background carrier wave

To better understand the dispersive properties of the nematic system given by the dispersion relation (3.6), we consider its long- and short-wave expansions. Expanding (3.6) in powers of *k*≪1 and retaining terms up to *O*(*k*^{5}) we have
*k*≪1, but that *νk*^{2}=*O*(1) or *νk*^{2}≪1, which generally does not hold true even for reasonably small wavenumbers *k* due to the very large value of the non-locality *ν*. Nevertheless, as we shall see, expansion (3.7) captures some key qualitative features of the full dispersion relation. Now looking at the short-wave asymptotics of (3.6), we obtain that for strong non-locality *ν*≫1
*ω*_{kk}<0, while for large wavenumbers *ω*_{kk}>0. Thus, the full dispersion relation (3.6) is non-convex, which has important physical consequences as it implies the possibility of resonance between long and short waves and hence the generation of short-wave radiation propagating ahead of the DSW. The effect of resonant radiation generation by solitary waves in equations with higher order dispersion is well known in the context of gravity–capillary waves (see [42,48–50] and references therein). There is also abundant literature on radiating solitary waves in nonlinear optics (see [52,43] and references therein). However, the counterpart of this for DSW theory is yet to be developed. A few existing notable contributions include numerical investigations of radiating DSWs described in the monograph [51] and the recent papers [12,14–17] on the effects of higher order dispersion on NLS DSWs in the context of nonlinear optics.

In figure 2, a comparison between the full dispersion relation (3.6) and the fifth-order Taylor expansion (3.7) is shown for the physically realistic non-locality *ν*=200 [47]. It can be seen that (3.7) is a good approximation to the full dispersion relation in the limit of low wavenumber, as expected. However, owing to the large factor in front of the *k*^{5} term in the approximate dispersion relation (3.7) the low wavenumber expansion rapidly deviates from the exact dispersion relation as *k* increases. Nevertheless, it qualitatively captures the key feature of the full dispersion relation, its non-convexity, so can be used for qualitative predictions of the effects of non-locality on the nematic DSW behaviour. It is further seen from figure 2 that the full phase velocity *ω*/*k* is not monotone and has a minimum, which is also qualitatively captured by the long-wave dispersion relation (3.7). The corresponding nonlinear equation with this linear dispersion relation, the fifth-order KdV equation, will be derived in the next section.

Let us now look at the opposite, dispersionless limit of the nematic system (3.2)–(3.4), which is described by the hyperbolic system of shallow water type [1]

## 4. Fifth-order Korteweg–de Vries equation

It has been shown that the nematic system (2.1) and (2.2) reduces to the KdV equation in the limit of small deviations from a background level [53,22]. However, the physically large value of the non-locality *ν* [47] and the linked resonant wavetrain have major effects on the asymptotic analysis, which were not considered in this previous work.

Indeed, assuming *ν*≫1, but *νk*^{2}=*O*(1), one has to retain the fifth-order terms in the dispersion relation expansion (3.7), which implies the necessity of keeping the fifth-order dispersion term in the corresponding asymptotic KdV equation. The asymptotic reduction of the nematic equations to the KdV equation in the limit of small deviations from a background level *u*_{0} will then be revisited, taking account of the large value of the non-locality *ν*. This asymptotic KdV equation will be derived from the hydrodynamic form of the nematic equations (3.2)–(3.4). Let us expand the hydrodynamic variables as
*ϵ*≪1 is a measure of the deviation from the background *u*_{0}, here *ϵ*^{2}=*u*_{2}−*u*_{1}. Also, *ξ*=*ϵ*(*x*−*Uz*) and *η*=*ϵ*^{3}*z* are the usual stretched variables used to derive the KdV equation [1]. We also assume that all corrections to the equilibrium state *ρ*=*ρ*_{0}, *v*=0, *θ*=*ρ*_{0}/*q* vanish as

Substituting expansions (4.1) and (4.3) into the director equation (2.2), we obtain at *O*(*ϵ*^{2})
*O*(*ϵ*^{4})
*νϵ*^{2}*θ*_{2ξξ}/2*q* is formally *O*(*ϵ*^{2}) and should appear at next order in the expression for *θ*_{3}, as in [53]. However, this implicitly assumes that *ν*=*O*(1), which is not the case for experimental values of *ν*. Hence, this term will be retained at *O*(*ϵ*^{4}). Treating *νϵ*^{2}*θ*_{2ξξ}/2*q* as a correction, equation (4.5) can be solved for *θ*_{2} to give
*P*_{2} can be of *O*(*ν*), making the last term *O*(*ν*^{2}*ϵ*^{2}).

Substituting expansions (4.1)–(4.3) into the ‘mass’ and ‘momentum’ equations (3.2) and (3.3), we have at *O*(*ϵ*^{3})
*θ*_{1}. Compatibility between these two equations for *V* _{1} and *P*_{1} then gives the coordinate velocity *U* as
*U*=*c* from the long-wave expansion (3.7) of the linear dispersion relation.

Similarly, at *O*(*ϵ*^{5}) the mass and momentum equations (3.2) and (3.3) give

It was shown in [22,53] that substituting the leading order part of (4.6) (the terms in brackets) into (4.10) and combining it with (4.7) and (4.9) leads to the KdV equation. We now need to extend this derivation by including the higher order terms of (4.6). The problem we encounter is with the computation of the last term in (4.6) as the correction *P*_{2} cannot be computed separately at order *O*(*ϵ*^{5}), leading to equations (4.9) and (4.10), and a higher order approximation is required. This difficulty can be circumvented by suggesting a suitable * ansatz* for

*P*

_{2}which is compatible with (4.9) and (4.10). Let

*α*is a constant. Then substituting (4.6) and (4.11) into (4.10) we obtain, on using (4.7),

Substituting (4.6), (4.11) and (4.12) into (4.9), we obtain the fifth-order KdV equation for *P*_{1}
*ξ* and *η* one has to replace (*ω*−*kc*)→*ϵ*^{3}*ω*, *k*→*ϵk* in (3.7) to make the comparison). We note that if the substitution (4.11) were not compatible with equations (4.9) and (4.10), it would not be possible to obtain agreement for both dispersive terms in (4.13) with the expansion of the nematic dispersion relation (3.7) using the single fitting parameter *α*.

The fifth-order KdV equation (4.13) differs from that found in [53,22] due to the *P*_{1ξξξξξ} term, which arises at this order as *ν* is large. The polarity of the solitary wave solution of the fifth-order KdV equation (4.13) depends on the sign of the coefficient of the *P*_{1ξξξ} term. It is then clear that in the non-local regime with *ν* large the solitary wave solution of the defocusing nematic equations (2.1) and (2.2) is a bright solitary wave, rising above a background level, rather than the usual dark solitary wave of the defocusing NLS equation, which the nematic equations become in the limit *ν*→0.

Although the fifth-order KdV equation (4.13) has a limited range of validity as an asymptotic, quantitative model for nematic DSWs, it provides major qualitative insight into their dynamics by capturing the effect of resonant radiation. To illustrate this, we solved numerically the normalized fifth-order KdV equation
*γ*>0. Equation (4.14) has been derived in several physical contexts, including magnetoacoustic waves and capillary–gravity waves of small amplitude when the Bond number is close to, but just less than,

Let us consider the fifth-order KdV equation (4.14) with the initial condition *w*=0, *x*>0 and *w*=*w*_{0}, *x*<0, so that a DSW is generated. Owing to the non-convexity of the dispersion relation for (4.14), there is the possibility of energy exchange between long and short waves propagating with the same phase velocity (figure 2), so this DSW is expected to generate a resonant linear wavetrain propagating ahead of it [51]. Such a radiating KdV DSW is displayed in figure 3. The solution shown in this figure has strong similarities to the radiating nematic DSW solution of figure 1*a*. However, the resonant wavetrain of the nematic solution is more uniform, which is due to the smoothing effect of the large non-locality *ν*.

In conclusion, we note that, although the known theory of radiating solitary waves provides some intuition as to the counterpart radiating DSW solution, the major contrasting feature of radiating DSWs is the fact that the lead solitary wave of the radiating DSW remains steady, while an isolated radiating solitary wave is intrinsically unsteady due to the radiation carrying away the solitary wave's energy [48].

## 5. Dam break problem for the nematic system

The solution of the Riemann problem (2.1)–(2.4) for the nematic system generically consists of three distinct parts: a rarefaction wave, a (bright) DSW and a radiative wavetrain (figure 1*a*). The rarefaction wave was analysed in [22], so below we only briefly outline the relevant results. Our main attention in this section will be on the DSW on the intermediate level *u*_{2} and the resonant wavetrain generated by it.

### (a) Rarefaction wave

The solution displayed in figure 1*a* shows that there is an expansion wave linking the initial level *u*_{3} behind the DSW to the level *u*_{2} on which the KdV DSW sits. This expansion wave solution has already been determined [22], so only the relevant details will be given here.

The expansion wave linking the initial level *u*_{3} behind to the intermediate level *s*_{+} is the velocity of the lead soliton of the KdV DSW, which lies at the leading edge of the intermediate shelf. The simple wave solution (5.1) linking the initial level *u*_{3} and the intermediate shelf *u*_{2} will be used for the comparisons with numerical solutions in §6.

The velocity *s*_{+} of the leading edge of the DSW on the shelf *u*_{2} needs to be determined. In contrast to the KdV and defocusing NLS equations, this velocity is not determined by the conservation of the Riemann invariant on *C*_{−} across the DSW [28], but by the classical shock jump condition. We note here that the occurrence of the classical shock conditions in a conservative dispersive hydrodynamics was observed earlier in numerical simulations of large amplitude shallow water DSWs [55] and optical DSWs in photorefractive media [56] (see also [57]) and, very recently, in the context of radiating dispersive shock waves governed by the defocusing NLS equation modified by the third-order dispersion [12,14,15,17]. This remarkable generic phenomenon requires further analytical study.

The non-dispersive equations (3.9)–(3.11) have the jump conditions [1]
*ρ*=*ρ*_{1} and *v*=0 and behind the shock, *ρ*=*ρ*_{2} and *v*=*v*_{2}. Eliminating between these equations gives
*v*_{2} in terms of the intermediate level *u*_{2}. Matching this and (5.4) then gives that this intermediate level *u*_{1}≤*u*_{2}≤*u*_{3}. For the particular case *u*_{1}=0, *u*_{1}→*u*_{3}, *u*_{2}→(*u*_{3}+*u*_{1})/2, which is the value obtained by conservation of the Riemann invariant (3.13) on *C*_{−} [22].

### (b) Dispersive shock wave: lead solitary wave

In previous work [22], the DSW solution of the KdV equation [25,26] was used for the DSW on the intermediate shelf *u*_{2}. While this was found to give good agreement with numerical solutions for values of *u*_{1} near *u*_{3}, significant disagreement was found for values of *u*_{1} away from *u*_{3}. As discussed above, this is due to the velocity *s*_{+} of the front of the full nematic DSW not being well determined by the velocity of the leading edge of the asymptotic KdV DSW. The reason for this behaviour is that the DSW is subject to radiative losses due to the resonance with the co-propagating linear short wavelength waves, resulting in a rapidly oscillating wavetrain shed ahead of the DSW. For small initial steps, the radiating DSW is described in the framework of the fifth-order KdV equation (4.13). However, for general jumps the full nematic system should be used due to the fifth-order KdV equation not being accurate in capturing large wavenumber dispersive behaviour (figure 2).

We now need to relate the shock velocity *s*_{+} to the amplitude of the lead soliton of the DSW. Since the solitary wave solution of the full nematic system is not available, as an approximation we shall use the soliton solution of the standard KdV equation, that is (4.13) without the fifth derivative term. On noting the scalings in expansions (4.1)–(4.3) and equating *s*_{+} given by (5.4) to the lead soliton velocity, this gives
*a* on identifying *u*_{0}=*u*_{1}.

### (c) Resonant wavetrain

Let us now consider the wavetrain ahead of the DSW. This wavetrain is generated due to a resonance between the long-wave oscillations in the DSW and co-propagating short wavelength waves, as implied by the non-convexity of the linear dispersion relation (3.6) and discussed in §4.

In determining the structure of the resonant wavetrain, we refer to figure 1*a* in which one can observe three regions of distinctly different behaviour: region (i) provides a transition from the lead soliton of the DSW to region (ii) which contains the (almost) uniform extended middle part of the wavetrain and the front region (iii) which brings the wavetrain down to the constant level *u*=*u*_{1} and *θ*=*Θ*_{1}, where *Θ*_{1}=|*u*_{1}|^{2}/*q*, see (2.4).

We start with the middle region (ii) for which the director *θ* is close to *Θ*_{1}, *θ*−*Θ*_{1}=*O*(*ν*^{−1}) and the wavenumber *k* is *O*(1). Hence, the asymptotic dispersion relation (3.8) applies and the dispersion relation for the resonant wavetrain is
*θ*=*Θ*_{1}:

We assume that the main resonance is with the lead soliton of the DSW, which we approximate by the KdV soliton (5.7). Matching the phase velocity to the lead KdV soliton velocity (5.4), we have
*c*_{g} [1], which is

The wavenumber (5.12) is real if *u*_{1}≤*u*_{1c}, where *u*_{1c} is the solution of
*u*_{1} above *u*_{1c}, there is only a transient wavetrain ahead of the DSW [1]. This existence of a critical *u*_{1} above which there is no resonant wavetrain is in agreement with previous work [22] in which the critical value was found to be *u*_{3}=1, *q*=2 and *ν*=200 numerical solutions give the critical value *u*_{1c}=0.69 [22]. For these parameter values, the new modulation value (5.14) *u*_{1c}=0.648 is slightly below the numerical cut-off, while the previous modulation value *u*_{1} range of about 0.1.

Above the critical value (5.14), the resonant wavetrain ceases to exist. The DSW on the intermediate level *u*_{2} then becomes the standard KdV-type DSW and the approximate solution of [22] holds. The amplitude and velocity of the lead soliton of the DSW are, then cf. (5.4) and (5.6),

As equation (5.10) is linear, it does not allow the determination of the resonant wavetrain amplitude. For that, one needs to go beyond the approximation *θ*=*Θ*_{1} in the wavetrain and include the (significant) variations of the director in the transition region (i) between the DSW and the uniform wavetrain region.

As the phase velocity of the resonant wavetrain is the same as the (classical shock) velocity *s*_{+} (5.4) of the lead soliton of the DSW, to determine the solution for the wavetrain in the transition region, we will use the moving coordinate *ζ*=*x*−*s*_{+}*z*. By inspection of the numerical solution of figure 1*a*, it is reasonable to assume that the approximate director solution in the transition region is given by the lead solitary wave of the DSW, so that from equations (4.4) and (5.7) of the KdV expansion of §§4 and 5b, we have
*σ*(*ζ*) and *u*_{r}(*ζ*,*z*) are the phase correction due to the variable coefficient and the wavetrain amplitude, respectively. To be consistent with the director (5.16) *u*_{r}=*O*(*ϵ*^{2}) as it is proportional to the jump height *u*_{2}−*u*_{1}. Substituting (5.16) and (5.17) into the electric field equation (2.1), we have
*σ*(*ζ*) so that the relation
*ϵ*, we obtain from (5.18) the equation for the variation of the wavetrain amplitude in the transition region as
*σ*′′ is higher order in *ϵ* (since

Using the numerical solution (figure 1*a*) and the soliton solution (5.7) as a guide to the structure of the transition region we shall look for the solution of equation (5.20) for *u*_{r} as fast (scaled as *O*(1)) oscillations with a slowly varying (scaled as *O*(*β*^{−1})) envelope. We then seek a WKB solution of the form
*X*=*βζ* and *Z*=*βz*. This WKB expansion is valid if *ν* is large. The first inequality is due to using the first two terms of the KdV expansion (4.3) for the director (5.16) and the second is required for the validity of the WKB form (5.21). Substituting the WKB form (5.21) into equation (5.20) gives the eikonal equation

We note that the group and phase velocity argument gave that as the resonant wavetrain approaches the wavefront at *x*=*c*_{g}*z*, it becomes a uniform wavetrain of wavenumber *k*_{r} and frequency *k*_{r}→*s*_{+}. This is expected as the group velocity of the front of the resonant wavetrain is *k*_{r}. When the velocity of the lead soliton of the KdV DSW is greater than the group velocity, the wavetrain cannot propagate away from the DSW. There is then no upstream resonant wavetrain, with only a small amplitude transient being present [22].

To solve the transport equation (5.23), the resonant wavetrain leading the KdV DSW must be matched to the intermediate shelf, so that *W*=*W*_{0}=*ϵ*^{2}=*u*_{2}−*u*_{1} at *X*=0 on noting the full solution (5.17) for *u*. Then using the eikonal equation solution (5.24), the solution of the transport equation (5.23) is
*x*=*c*_{g}*z* is approached, so that the total height of the envelope of the resonant wavetrain in the region (ii) is given by

Finally, we describe region (iii) of the resonant wavetrain which brings it down to the initial level *u*_{1} (figure 1*a*). In the region of this front, as for the uniform middle region (ii), we approximate *θ* by *θ*=*Θ*_{1}=|*u*_{1}|^{2}/*q*, so that the linear equation (5.10) holds. If we use a moving coordinate *ζ*_{g}=*x*−*c*_{g}*z* moving with the velocity of the front, the electric field is governed by
*u*_{f} is the solution of
*u*_{f}|=*a*_{r}−*u*_{1} at *ζ*_{g}=0. The linear equation (5.30) can be solved using Laplace transforms to give the Fresnel integral solution

## 6. Comparison with numerical solutions

In this section, full numerical solutions of the nematic equations (2.1) and (2.2) will be compared with the modulation theory solutions of §5a–c. The numerical solution of the electric field equation (2.1), which is of NLS-type, was obtained using the pseudo-spectral method of Fornberg & Whitham [26], modified to improve its accuracy and stability [58], but without the boundary damper due to the non-zero boundary conditions. These improvements include using a fourth-order Runge–Kutta scheme to propagate forward in *z*, resulting in higher accuracy, in Fourier space, rather than in real space, resulting in improved stability [58]. The numerical solution of the linear director equation (2.2) was obtained using a spectral method [59]. This numerical scheme is discussed in [60]. For the numerical solutions of this work, 32 768 points were used for the FFT with a *x* domain of length 8192.0 and a *z* step of *Δz*=0.002. The *x* domain was chosen long enough so that the waves at the numerical boundaries generated by periodicity were far from the region of interest. Finally, the initial condition (2.3) was smoothed using *x* derivatives, with *W*=1 found to be suitable.

Figure 4 shows a comparison between the numerical solution of the nematic equations (2.1) and (2.2) for *u*_{3}=1.0 and *u*_{1}=0.5 at *z*=300 for *q*=2 and *ν*=200. For clarity, in these figures only the upper envelope of the resonant wavetrain (5.17), (5.21) and (5.25) and the upper envelope of the Fresnel front (5.31) are shown. It can be seen that there is very good agreement in general between the numerical solution for the electric field |*u*| and the modulation theory solution of §5a–c. In particular, there is excellent agreement for the position of the lead soliton of the DSW, which is the same as that of the trailing edge of the resonant wavetrain. This is in contrast to the result of previous work [22] in which this position was determined by the velocity *x*=300 for the parameters of figure 4, noting that the numerical position is *x*=247.5. It is then clear that the shock velocity (5.4) determined from the shock jump conditions for the non-dispersive ‘shallow water’ equations (3.9)–(3.11) and giving *x*=265.9 for the lead soliton at *z*=300 yields much better agreement with the numerical solution for the position of the leading edge of the DSW than the velocity determined by the KdV DSW solution. The differing length scales of the KdV DSW (*O*(1)) can be clearly seen. The major disagreement is that the front of the numerical resonant wavetrain has more structure than the linear Fresnel integral solution of §5c. However, the Fresnel integral solution gives the correct spatial extent of the transient front of the resonant wavetrain. Furthermore, if the Fresnel integral solution is shifted so that it starts ahead of the rise in the numerical front, it is in very good agreement with the numerical front.

The other notable disagreement between the numerical and analytical solutions is the amplitude of the lead soliton of the DSW. The amplitude of the DSW in the electric field *u* is generally under-predicted by the KdV theory of §5b, which was based on the classical shock speed (5.4). However, this approximation yields good agreement for the DSW in the director *θ*, given in the KdV approximation by equation (5.16). This is in contrast to the results of [22] for which the standard DSW solution of the KdV equation was used to determine the DSW on the intermediate level *u*_{2}. The results of [22] strongly over-predicted the height of the bore in the director, this major discrepancy being fixed in this theory. Finally, it can be seen that under the resonant wavetrain there is a slight rise in the director above *θ*=*Θ*_{1} due to *O*(*ν*^{−1}) corrections in the asymptotic expansions. These higher order corrections will be dealt with in future work based on a full description of a resonantly radiating DSW.

The agreement between the modulation theory and numerical solutions is further quantified in figure 5. Figure 5*a* shows a comparison of the height (background plus amplitude) of the lead soliton of the DSW as given by numerical solutions and the modulation solution (*a*_{s}=*ϵ*^{2}*A*+*u*_{1}), using (5.6) for the amplitude below the cut-off (5.14) and (5.15). The choice of the total height rather than amplitude for the comparisons is due to the soliton background being not clearly defined in the numerical solutions (figure 1*a*). Furthermore, the amplitude of the lead wave oscillates slightly due to its interaction with the resonant wavetrain, so the figure shows the average amplitude. The numerical solution clearly shows the predicted different DSW behaviours above and below the resonant wave cut-off, which was not predicted in [22] for which the height was the constant value (5.15) for the whole range of *u*_{1}. Thus, one can see that the KdV soliton height based on the classical shock wave speed is in broad agreement with the numerical values. The appropriateness of using the classical shock wave velocity to determine the intermediate shelf height (5.5) is quantified in figure 5*b*. It can be seen that (5.5) is in excellent agreement with the numerical height, except for a slight discrepancy as *u*_{1}→0. This is due to the intermediate shelf disappearing as the dam break solution for *u*_{1}=0 is approached [22].

Figure 5*c* shows a comparison between the height of the resonant wavetrain obtained numerically and the modulation solution height (5.27). There is excellent agreement between these heights, except towards the cut-off near *u*_{1}=0.7. This is due to the discrepancy between the numerically found cut-off and the modulation theory prediction.

Finally, figure 5*d* shows a comparison for the leading and trailing edges of the resonant wavetrain. It can be seen that there is excellent agreement for the position of the trailing edge, even up to the cut-off. Previous work [22] predicted the constant (i.e. independent of *u*_{1}) velocity (5.15) for the trailing edge which was defined by the value *u*_{3} alone. Figure 5*d* clearly shows that the present theory based on the classical shock speed for the leading edge of the DSW is in much better agreement. The agreement for the leading edge is reasonable above *u*_{1}=0.5 as the cut-off is approached, but is poor as *u*_{1} decreases. There are a number of reasons for this. The position of the trailing edge is clearly defined by the peak of the lead soliton of the KdV DSW. While the theory of §5c predicts a precise location for the leading edge of the resonant wavetrain, it can be seen from figures 1*a* and 4 that this is not the case for the numerical solution. There is no clean boundary between the resonant wavetrain and its front. There is an extended transition between the two. For the comparison of figure 5*d*, the start of the hump before the monotonic decrease of the front was chosen as the leading edge position of the wavetrain. Finally, the present modulation solution underpredicts the cut-off point for the resonant wavetrain (at *u*_{1}=0.648 compared with the numerical value 0.69 for the parameter values of figure 5), leading to the disagreement as the cut-off is approached.

Part of the reasons for these discrepancies is due to the major simplifying assumption adopted in the present theory, according to which the generation of the wavetrain is dominated by a single resonance with the lead soliton of the DSW. A more advanced theory including internal resonances with the other components of the DSW is needed to achieve better agreement with numerical solutions.

There is one more feature of the resonant wavetrain which complicates its analysis for large initial jumps. As the initial level ahead *u*_{1} decreases the electric field *u* eventually vanishes at a point, termed the vacuum point [37]. For sufficiently large initial jumps, the vacuum point occurs within the resonant wavetrain, so that the lower envelope becomes non-monotone. It was shown in [37] that for the defocusing NLS DSW there is a singularity in the phase *v* at the vacuum point itself. Although the resonant wavetrain for the nematic system (2.1) and (2.2) is asymptotically described by the linear equation (5.10) rather than the defocusing NLS equation, numerical simulations show that the vacuum point in the wavetrain has qualitatively similar properties to the vacuum point arising in the large amplitude NLS DSW [37]. In particular, such a DSW has a non-monotone lower envelope (figure 6*a*) and exhibits a phase singularity at the vacuum point (figure 6*b*). The WKB solution of §5c gives that the lower envelope of the resonant wavetrain has height
*u*_{3}=1.0, *ν*=200 and *q*=2, it is found that *a*_{l} vanishes when *u*_{1}=0.2416. Numerical solutions of the nematic equations (2.1) and (2.2) show that for these parameter values a vacuum point first occurs when *u*_{1}=0.22. A full analysis of the solution after the vacuum point is reached is beyond the scope of this paper. Full Whitham modulation equations would be required for a proper analysis after the vacuum point [37].

## 7. Conclusion

The Riemann problem for the equations governing the propagation of a coherent optical beam in a defocusing nematic liquid crystal has been studied. It was found that in the highly non-local limit the DSW, which comprises a major part of the Riemann problem solution, is drastically different to the DSW solution of the defocusing NLS equation, to which the nematic equations reduce in the small non-locality limit, that is *ν*→0. There are two major differences: (i) the nematic DSW is of positive polarity with a bright soliton at its leading edge; (ii) it is preceded by a short wavelength resonant wavetrain. To clarify this structure, it was shown that in the limit of small deviations from a background, the nematic equations reduce to a KdV equation with a fifth-order derivative, the Kawahara equation. This fifth-order KdV equation is known to have a resonance between its solitary wave solution and linear radiation. The present work shows that this resonance extends to a resonance between the DSW and linear radiation. A modulation theory was developed to derive solutions for the resonant wavetrain and its front. In contrast to previous work [22], it was found that the leading edge of the DSW was determined by the classical shock jump condition, which is non-standard for DSWs [28]. Excellent agreement was found between the major part of the modulation theory solution and full numerical solutions of the nematic equations. However, there are some discrepancies. Part of the observed discrepancies can be addressed by applying a more complete modulation theory for the DSW description. When the higher order dispersion term is a small perturbation, as in the Kawahara equation (4.14) with *γ*≪1, Whitham theory for perturbed integrable equations [61,62] provides an appropriate analytical framework for the description of the DSW evolution.

The present work leaves open a number of issues. As already discussed, resonances between DSWs and radiation is an issue which has received little attention to date with the theory and solution methods only starting to be developed [51,12], in contrast to the corresponding resonant interaction between solitary waves and radiation (e.g. [43,48–50,52]). As the nematic equations are generic and similar equations arise in a number of fields, this resonant interaction deserves in-depth treatment. A proper analysis of possible resonances between DSWs and radiation is an open question which deserves further treatment. This will be the subject of further work.

## Ethics

This work does not involve any experiments on humans or animals.

## Data accessibility

This paper does not use any experimental data.

## Author' contributions

G.A.E. and N.F.S. contributed equally to the analysis of this paper, with N.F.S. performing the numerical calculations.

## Competing interests

The authors have no competing interests.

## Funding

G.A.E. and N.F.S. were supported by the London Mathematical Society under the Research in Pairs Grant 41421. G.A.E. was also partially supported by the Royal Society International Exchanges Scheme Grant IE131353.

## Acknowledgements

G.A.E. thanks Mark Hoefer for discussions concerning the Kawahara equation.

- Received September 7, 2015.
- Accepted February 4, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.