arXiv:0704.1212v1  [physics.optics]  10 Apr 2007
CSHOCK-II Dr.Paul.Kinsler@physics.org
Optical carrier wave shocking: detection and dispersion
P. Kinsler, S.B.P. Radnor, J.C.A. Tyrrell, and G.H.C. New
Blackett Laboratory, Imperial College London, Prince Consort Road, London SW7 2BW, United Kingdom.
(Dated: November 20, 2018)
Carrier wave shocking is studied using the Pseudo-Spectral Spatial Domain (PSSD) technique.
We describe the shock detection diagnostics necessary for this numerical study, and verify them
against theoretical shocking predictions for the dispersionless case. These predictions show Carrier
Envelope Phase (CEP) and pulse bandwidth sensitivity in the single-cycle regime. The flexible
dispersion management offered by PSSD enables us to independently control the linear and nonlinear
dispersion. Customized dispersion profiles allow us to analyze the development of both carrier self-
steepening and shocks. The results exhibit a marked asymmetry between normal and anomalous
dispersion, both in the limits of the shocking regime and in the (near) shocked pulse waveforms.
Combining these insights, we offer some suggestions on how carrier shocking (or at least extreme
self-steepening) might be realised experimentally.
I. INTRODUCTION
The self-steepening of an optical pulse envelope was
first studied by DeMartini et al. in 1967 [1], and is a
well-known phenomenon associated with self-phase mod-
ulation (SPM). Surprisingly however, the possibility of
self-steepening of the optical carrier wave1 was considered
even earlier in a 1965 paper by Rosen [2], who showed
that, for a third order χ(3) nonlinearity (and under suit-
able conditions), a shock (or field discontinuity) develops
in a finite distance. This latter phenomenon received lit-
tle attention for more than 30 years, until it was revisited
in the 1990s by Moloney et al. [3, 4], who performed Fi-
nite Difference Time Domain (FDTD) simulations of the
process. Three dimensional FDTD simulations of carrier
shocking have recently been performed by Trillo et al.
In the present paper, we investigate carrier wave shock-
ing in χ(3) nonlinear materials where SPM is accompa-
nied by the generation of (odd) higher harmonics. For
the dispersionless case we have generalized earlier pre-
dictions based on the method of characteristics (MOC)
to allow for arbitrary initial waveforms. This allows us,
for example, to predict the carrier envelope phase (CEP)
sensitivity of the shocking distance for optical pulses, as
well as the dependence on pulse length.
Our primary interest in this paper is the effect of linear
dispersion on carrier shock formation. Since no analytic
solutions exist in this case, we are forced to rely on nu-
merical simulations. In any discussion of carrier wave
shocking, it is important to distinguish between three
related concepts: the physical system, the mathemati-
cal model, and the numerical model. The mathematical
model is an approximation to the physical system, while
∗Electronic address: Dr.Paul.Kinsler@physics.org
1 Note that in this paper we use the term “carrier” to denote
all the oscillations of the field, and do not use its other sense,
i.e. that of fixed-frequency oscillations as used in a envelope and
carrier representation of a pulse.
the numerical model is an approximation to the mathe-
matical model. Actual discontinuities occur only in the
mathematical model, and it is these to which the idea of
shock formation refers. In the physical system, disconti-
nuities are prevented by phenomena not included in the
mathematics, while numerical codes inevitably fail as a
mathematical discontinuity is approached. The indicator
of an imminent “shock” is the rapidly increasing gradi-
ent (steepening) of the optical carrier. In the numerical
model, we see numerical symptoms generated by extreme
self-steepening, and these correlate with the onset of the
discontinuities in the mathematical model. Under these
circumstances, it has been necessary for us to develop
a quantitative numerical test of “shock formation”. We
have found the most satisfactory diagnostic to be “Local
Discontinuity Detection” (LDD), which gives results that
are in good agreement with theoretical MOC predictions
based on the mathematical model. LDD provides a clear
numerical measure of the rapid steepening that precedes
the appearance of a discontinuity in the mathematical
representation. We continue to use the term “shock” to
describe the situation where a mathematical shock is im-
minent.
In our simulations, we exploit the flexibility in disper-
sion management offered by the Pseudo-Spectral Spatial
Domain (PSSD) technique [7] to study carrier shock for-
mation for a range of simple dispersion profiles, and de-
termine the degree of phase mismatch that carrier shock-
ing can tolerate. As expected, we find that shocking oc-
curs when the nonlinearity dominates the linear disper-
sion. However, it emerges that the process is asymmet-
ric, with anomalous dispersion being far more conducive
to shocking than normal dispersion. Hence, the fact that
the LDD scheme does not detect carrier shocks in simula-
tions involving (normally dispersive) fused silica, even at
powers equal to its damage threshold, is neither surpris-
ing nor necessarily discouraging. After all, anomalously
dispersive materials could potentially be engineered. Fur-
ther, our results also relate to how one might perform
carrier shaping (as opposed to carrier steepening), a pro-
cess that has some interesting applications.
1 Kinsler-RTN-2007
http://arxiv.org/abs/0704.1212v1
mailto:Dr.Paul.Kinsler@physics.org
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
After briefly describing our simulation methods in
section II, we consider dispersionless shocking and the
method of characteristics in III, and our LDD shock de-
tection scheme in IV. Then we discuss the effect of disper-
sion on carrier shocking in V, followed by our numerical
results in VI. In VII we consider the potential relevance
to experimental detection of carrier steepening and/or
shocks. Finally, in VIII, we present our conclusions.
II. SIMULATION METHODS
The PSSD method [7, 8] offers significant advantages
over the traditional FDTD and Pseudospectral Time-
Domain (PSTD) [9] techniques for modeling the propa-
gation and interaction of few-cycle pulses. Run times are
generally faster, and PSSD also offers far greater flexi-
bility in the handling of dispersion. Whereas FDTD and
PSTD [9] propagate fields E(z), H(z) forward in time,
PSSD propagates fields E(t), H(t) forward in space. It
is important to keep this difference in mind when com-
paring our results to those in [9]. Under PSSD, the en-
tire time-history (and therefore frequency content) of the
pulse is known at any point in space, so arbitrary disper-
sion incurs no extra computational penalty. In contrast,
the FDTD or PSTD approaches use convolutions incor-
porating time-response models for dispersion.
We apply the PSSD algorithm to two representations
of the field and source-free Maxwell’s equations in non-
magnetic media; the first uses the E and H fields, and
the second the directional fields ~G±(t) = αr(t)∗ ~Ex(t)±
βr ~Hy(t) [10]. Here the αr, βr include the (linear) per-
mittivity and permeability of the material (i.e. ǫ(t), µ).
These G± fields enable us to rewrite Maxwell’s equations,
and efficiently separate out the relevant forward-going
part of the field.
For an instantaneous χ(3) nonlinearity, the equations
for E and H in the 1D (plane wave) limit are
dHy(t; z)
ǫ0ǫr(t) ∗ Ex(t; z) + χ
(3)Ex(t; z)
dEx(t; z)
[µ0Hy(t; z)] . (2)
The G± field simulations usually assume G− = 0, and as
a result contain only forward traveling components. The
forward-only wave equation for G+ is
dG+(t; z)
βrαr(t) ∗G
+(t; z) + βrχ
(3)E(t; z)3
where it is most straightforward to calculate the non-
linear term by reconstructing E(t; z) from G+(t; z) in
the frequency domain using Ẽ(ω; z) = G̃+(ω; z)/2α̃r(ω),
since G− = 0. Notice the similarity between eqns. (1)
and (3), but that eqn. (3) propagates the field in a single
first order equation, rather than two.
Typical array sizes used in pulse simulations were
N = 214 covering a time window T = 200fs, and (spatial)
propagation steps were dz = 0.4cT/N ≈ 0.9nm. We en-
sured the stability of our integration using Orszag’s 2/3
rule [11], which involves setting the upper part of the
spectral range to zero. It is worth noting that changing
this cut-off, either by adjusting its position, or by using
a smoothed (rather than step-like) filter, made little dif-
ference to test simulations. The pulse profile used as an
initial condition was
E(t) = E0 sin(ω1t+ φ) sech(0.28ω1t/τ), (4)
where our standard parameters were ω1 = 2.356 ×
1015rad/s (i.e. λ = 800nm) with τ = 0.93̇. Such pulses
are rather short (since the number of cycles inside the
intensity FWHM is τ), but in fact the shocking distance
is only weakly dependent on the pulse width, with signif-
icant variation only appearing for pulses of a few cycles
or less.
We also performed CW simulations, for which we mod-
eled just a single cycle of the carrier. assisted by the pe-
riodic nature of the discrete Fourier transform. For these
we used array sizes of N = 210, and the time window was
set by the period of the field oscillations.
Our default value of nonlinear strength was χ(3)E20 =
0.02, which is comparable to that in fused silica at an in-
tensity of 0.7×1014 W/cm2; our χ(3) parameter is equiv-
alent to η in Rosen [2], and a in Gilles et al.[4]. We use an
instantaneous nonlinearity, since our primary interest is
linear dispersion, and that is a far more significant effect
than the nonlinear time response.
III. CARRIER WAVE SHOCKING
As an introduction to the process of carrier wave shock-
ing, fig. 1 shows the profile of a pulse propagating in a
dispersionless χ(3) medium, just before a shock occurs.
The nonlinearity gives rise to a nonlinear index of refrac-
tion n2E
2, the effect of which is to increase the effective
refractive index in the more intense regions of the pro-
file. This reduces the phase velocity at the peak of each
oscillation with respect to the rest of the waveform, and
causes the slope on the trailing edges to increase dramat-
ically.
The effects seen in fig. 1 are associated with the gen-
eration of third and higher harmonics, although the har-
monic components are not particularly strong even when
a shock is about to occur. As an example, fig. 2 shows
how the harmonics build up as shocking is approached.
The ω4 scaling of the intensity spectrum in the figure ex-
aggerates the contribution of the higher orders, and has
been chosen for illustrative purposes. Notice that the
profile becomes nearly flat (i.e. the spectrum falls off as
the 4th power of the frequency) just before shocking is
registered at around 4.3µm.
If we write
E = A(t) [sin(ω1t) + γ cos(3ω1t+ ψ)] , (5)
2 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
5 10 15 20
time(fs)
FIG. 1: The profile of a few-cycle optical pulse just prior to
shocking in the dispersionless limit. The larger oscillations
in the centre of the pulse undergo more self-steepening than
those in the wings. The standard pulse parameters were used.
FIG. 2: Development of the heights of the scaled harmonic
peaks as a pulse approaches the (LDD) shocking distance of
4.3µm in the dispersionless case. Each line perpendicular to
the ω/ω1 axis corresponds to the (scaled) contribution from
that spectral peak. Note that the viewpoint has been rotated
so that the contribution from the fundamental is to the right.
The initial pulse contained about 33 cycles (τ = 33).
we find that an appropriate choice of A(t), with γ =
0.1 and ψ = 0, gives us a passable match to fig. 1.
Note that this choice of ψ corresponds to the phase of
third harmonic generation (THG) under index matched
conditions, i.e. where n0 = n(ω1) = n(ω3).
Rosen’s original paper [2] used the MOC to predict
the formation of a value discontinuity in the field at cer-
tain points within the profile. If the displacement of the
dispersionless medium is written
D = ǫ0
E + χ(1)E + χ(3)E3
, (6)
he showed that the wave equation for E is
1 + χ(1)
) ∂2E
+ χ(3)
, (7)
and that the associated equation governing the charac-
teristic lines of E is
+ v(E)
= 0. (8)
Here, the velocity v(E) is given by
v(E) =
(ǫr + 3χ(3)E2)1/2
, (9)
where ǫr = 1 + χ
(1) = n0
2 is the (relative) dielectric
constant and n0 the linear refractive index.
v − dv
FIG. 3: Method of Characteristics. Two points A and B
on the field profile, separated initially by a time difference dt
travel at different speeds v − dv and v, and meet at point C.
Using eqn. (9) along with the construction shown in
fig. 3, we can derive a simple formula for the distance to
shocking. The figure shows two characteristics AC and
BC, originating from points A and B, and converging
towards a shock at C after a distance of L. The intensity
at A is higher than at B, so the speed associated with
AC (represented by its gradient) is lower than that of
BC. This means that at C, the field has two values, and
a discontinuity has formed. From the geometry of the
figure, it is easy to show that
(c/n0)
, (10)
where t, v = c/n0, and L = vt are respectively time,
speed and distance.
On the other hand, differentiating eqn. (9) leads to
3cχ(3)
n2o + 3χ
(3)E2
, (11)
and this combined with eqn. (10) yields
1 + 3χ(3)E2/n2o
3χ(3)(−dE2/dt)
1 + 8n2E2/no
(−dE2/dt)
, (13)
where n2 = 3χ
(3)/8no is the material parameter deter-
mining the intensity induced refractive index shift as
2. For a given profile, a shock will occur first at
the point where −dE2/dt reaches its negative extremum.
3 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
We can therefore define the shocking distance as
Lshock =
1 + 8n2E2/no
(−dE2/dt)
(−dE2/dt)
for 8n2E
2/no ≪ 1.(15)
This formula is more general than that of either Rosen
[2] or Gilles et al. [4]. Notice in particular that the
parameter that controls the shock behaviour is not the
gradient of the field (dE/dt), but that of the field squared
(dE2/dt).
For a sinusoidal initial waveform E = E0 sin(ω1t+ φ),
it is easy to show that shocks form at
ω1t = −π/4 + jπ − φ, (16)
where j is an integer, and that the shocking distance is
Lshock =
3ω1χ(3)E
4ω1n2E
. (17)
The approach is readily extended to pulsed waveforms
of the type defined in eqn. (4). Analytical results can
be derived in the new situation in which eqn. (16) be-
comes a transcendental equation. However, the results
are cumbersome and it is simpler to scan the profile nu-
merically to determine the shocking parameters. It turns
out that the shocking distance for short pulses exhibits
an interesting sensitivity to the carrier envelope phase φ.
Whilst in the CW case all locations defined by eqn. (16)
were equivalent, for pulses, the one nearest the peak of
the envelope has a shorter Lshock than the others.
A set of results is displayed in fig. 4 where the shock-
ing distance is plotted as a function of the carrier phase
φ for sech profiles with different pulse widths τ . The
dotted line is for the case of a very broad pulse for which
Lshock is given by eqn. (17). The sharp peaks mark a
curve crossings where the shock location switches from
one point on the sinusoid to another. Notice that the
range of shocking distances increases for shorter pulses,
and that, unsurprisingly, the lowest values occur when φ
is around π/2 i.e. when the pulse has a cosine form.
One very important point to note from fig. 4 is that for
pulses containing more than a few cycles, the dependence
of the shocking distance on pulse width is very weak. We
will see this message repeated later in section V, with
shock regions being similar for both single cycle (τ = 1)
pulses and CW (τ = ∞) fields.
IV. SHOCK DETECTION
Since optical shock formation is directly associated
with regions of increasingly steep field gradient, any nu-
merical scheme, however sophisticated, is bound to fail at
some point in the process. We therefore want to recognize
when a shock is imminent, not only to avoid numerical
FIG. 4: MOC shocking distances as a function of phase,
for pulses as in eqn. (4), allowing for different pulse lengths.
We include the 1/2 cycle τ = 0.56 results to emphasize the
trend; but even for a long pulse (τ = 28), a small peak can
be seen just below φ/π = 0.01. The peak position is weakly
τ dependent. The + signs denote LDD shocking distances
obtained from simulations of the τ = 0.93̇ case.
problems, but as a means to estimate the distance at
which a discontinuity would occur in the mathematical
model.
One obvious symptom of impending numerical failure
is loss of energy conservation [7]. However, this does not
give an indication at the first instance of a shock form-
ing, but rather signals the accumulated effect of multiple
small numerical failures from many shocked regions.
A more physical strategy is to search for regions where
the field gradient dE/dt is large and increasing rapidly,
and to use this to predict the shocking point. A useful
variant, suggested by the MOC calculation in the previ-
ous section, is to use the value of −dE2/dt instead of the
gradient.
Overall, we find that the best method is Local Discon-
tinuity Detection (LDD), which is similar to techniques
used in other fields (see e.g. [12]). As the shock regime is
approached, narrow shoulders with associated points of
inflection appear within the regions of rapidly increasing
gradient. The procedure is therefore to scan the field pro-
file for the maximum gradient (of either E or E2) and, if
it occurs near a point of inflection, an incipient shock is
registered. An example of a pulse that has just triggered
the LDD diagnostic is shown in fig 5.
The LDD method requires two parameters. The first
determines the time scale used in determining whether a
point of inflection exists. For this, we pick the scale set
by our temporal grid and insist that the field gradients
calculated at three adjacent grid points have opposing
signs: either up-down-up, or down-up-down. The sec-
ond determines the maximum range allowed between the
maximum gradient and the point of inflection, and our
default value for this was 10 grid points. In our simula-
tions, we see that the position of the first detected shock
depends only weakly on this range. We can easily mini-
mize the small sensitivity to these parameters by holding
4 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
them fixed throughout any given set of simulations. As a
result, we have found LDD to be a sensitive and reliable
method of shock detection.
FIG. 5: The profile of a few-cycle optical pulse at the point
of shocking. The larger oscillations in the centre of the pulse
undergo more self-steepening than those in the wings. (a) The
whole pulse, with the LDD carrier shock circled (◦). (b) An
expanded view of the shock region, with its point of inflection
(at ∼ 9.3fs) being very close to the steepest gradient (at ∼
9.1fs). The standard pulse parameters were used.
Eqn. (17) predicts that the shocking distance should
increase linearly with the refractive index n0. We use
this to test the LDD diagnostic in fig. 6, where the ana-
lytical formula is compared with the results of numerical
simulations. This figure shows close agreement between
prediction (using eqn. (15)) and simulation, where the
approximation causes the MOC prediction to be reduced
by less than 0.1µm at n0 = 1. We see similar agreement
between simulations using the LDD method and the car-
rier phase sensitivity shown in fig. 4. The presence of
small systematic differences (as on e.g. fig. 4, 6) can be
easily understood, since the LDD diagnostic is (strictly
speaking) a test of the numerics, and is not a direct test
for the presence of a physical shock or mathematical dis-
continuity.
V. THE EFFECT OF DISPERSION ON
SHOCKING
The primary purpose of this paper is to understand
the principles of carrier shock formation in the presence
of dispersion. Some simple ideas about the role of dis-
persion can be understood from eqn.(5) using the insight
FIG. 6: The LDD shocking distance as a function of refractive
index, comparing the approximate MOC prediction from eqn.
(15) to PSSD simulations (+) for pulses with τ = 0.93̇. Other
PSSD simulation results from independent codes give very
similar results to those shown on the graph.
from eqns.(10) – (15), i.e. that the rate of change of E2
is the critical factor in shock development, rather than
that of E. If only the fundamental and third harmonic
components are considered (as in eqn.(5)), the effective
refractive indices are
∆nNL:1 = n1 + n2
1 + 3γ2
I1, (18)
∆nNL:3 = n3 + n2
3 + γ2
I1. (19)
Evidently, the relative phase velocity of the two waves is
affected by both linear and nonlinear dispersion, so the
phase ψ in eqn.(5) will vary accordingly as the pulse prop-
agates. Fig. 7, which shows how dE2/dt varies with ψ,
suggests that shocking is likely to be exacerbated when ψ
is small and positive, but moderated when ψ is negative.
Broadly speaking, the former case will be promoted by
anomalous dispersion and the latter by normal disper-
sion; however, the process is clearly complicated, since it
involves time-dependent phase shifts between the waves
and the interplay of linear and nonlinear dispersion. Of
course, since carrier shocking relies on the establishment
and maintenance of specific phase relationships between
a set of harmonics, it must be expected that strong dis-
persion of either sign will disrupt the shock formation
process. On the other hand, the simple argument that
has been offered suggests that shocking may be tolerant
to a degree of anomalous dispersion, but not to a simi-
lar amount of normal dispersion. In general, therefore, a
graph of shocking signature versus refractive index mis-
match might be expected to exhibit a shock region where
nonlinearity dominates dispersion (displaced in the direc-
tion of anomalous dispersion), surrounded by a shock-
free region where dispersion dominates the nonlinearity.
Moreover, if a dominant coherence length LC can be de-
fined, it is reasonable to expect shocking to occur when
this exceeds the characteristic SPM length (LSPM ). As
we shall see, all these features are borne out by the nu-
merical results.
5 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
FIG. 7: The effect of different time-lags between fundamental
and third harmonic on the steepness of the pulse. The field
profile is as in eqn. (5), with A(t) = A0 and γ = 1/16. (a)
Scaled maximum values of Mj = Max(dE
j/dt) as a function
of third harmonic phase offset ψ. We can see that the maxima
of dE/dt and dE2/dt occur at different ψ. (b) Position on the
pulse θj = ω1tj of the maxima plotted in the top frame. Note
that the kinks in (a) and (b) occur for offsets that give the
longest shocking distances, not the shortest.
In our numerical simulations, we make extensive use
of model dispersion curves, which enable refractive in-
dex differences between harmonics to be freely controlled,
and lead to uncomplicated boundaries between the shock
and no-shock regimes. The dispersion profiles shown
in fig. 8 duly contain either refractive index steps ∆n
at the midpoints between successive odd harmonics, or
smooth gradients δ that provide a similar net change.
In the simplest option, there is a single step at 2ω1,
in which case the dominant coherence length is clearly
LC = π/ |k3 − 3k1| = πc/3ω1 |n3 − n1|; we have also
tried a single step at 4ω1, which has a shorter LC . In
all cases, when n increases with frequency, the situation
corresponds broadly to normal dispersion, while decreas-
ing n corresponds to anomalous dispersion. The step
size (or gradient) is chosen to be comparable to values in
fused silica at ω1 = 1.5rad/fs, where the refractive index
differences between the lower harmonics are ∆n1,3 ≈ 0.06
and ∆n3,5 ≈ 0.12[4]. While we consider the case of fused
silica itself in Section VII, the results based on the model
dispersion characteristics are invaluable for understand-
2   n∆
ω1 ω3 ω5 ω7
FIG. 8: Types of refractive index profile used. The solid lines
show a single refractive index step midway between ω1 and
ω3, so that the fundamental is always phase mismatched from
its higher harmonics, although they remain perfectly phase
matched with each other. The dashed lines show multiple re-
fractive index steps, each of the same size, and always midway
between subsequent harmonics. The dotted lines show a lin-
ear refractive index gradient which gives the same mismatch
between subsequent harmonics as the multi-stepped case. We
do not show lines for the case of a single step at 4ω1 to avoid
cluttering the figure.
ing the essential principles of carrier shock formation in
dispersive media.
VI. RESULTS AND DISCUSSION
We will now analyze our numerical simulations of car-
rier shocking in the presence of dispersion on the basis of
the principles discussed above. Results are included for
both a CW wave and for a single-cycle pulse. The CW
case gives slightly wider shocking regions, but the differ-
ences are minor. This is because the effect of the pulse
envelope on the field amplitudes and gradients of the cen-
tral carrier oscillation is small, except when considering
sub-cycle pulses. In the results we present, the energy
conservation and LDD measures reveal the presence of
sharp boundaries between the shocking and non-shocking
regions.
Results for single refractive index steps at 2ω1 and 4ω1
are presented in figs. 9 and 10 respectively, whereas fig.
11 has a step midway between all harmonics. In both
single stepped cases, a useful coherence length can be
defined on the basis of the mismatch between the fun-
damental and the harmonic just above the step. In the
multi-stepped case, there is no easy way to define a dom-
inant coherence length.
Fig. 9(a) shows the pulse profiles in the anomalous
(negative) and normal (positive) single step cases, for the
smallest step size at which shocking did not occur. The
first obvious difference between the profiles is the oppo-
site walk-off direction of the third harmonic, although the
third harmonic contribution is hard to see in the anoma-
lous case. The second is that the pulse profiles exhibit
6 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
−20 −15 −10 −5 0
time(fs)
−0.1 −0.05 0 0.05 0.1
∆ n1,3
−0.1 −0.05 0 0.05 0.1
(b) X
FIG. 9: Carrier shocking for CW and pulsed (τ = 0.93̇) cases
with a single refractive index step at 2ω1. Frame (a) compares
chirped pulses at a 40µm propagation distance from the non-
shock region immediately outside of the shocked region; the
upper curve (X) is for the negative step, the lower (Y) for
the positive step. Frame (b) shows the correlation in the
region where LC . LSPM between energy conservation failure
(logarithmic right hand scale, dots) and the LDD detected
shocking distance (left hand scale, solid line). The dotted
line shows the LDD results for a CW field (τ = ∞).
distinctly different characteristics according to the sign of
the dispersion. Narrow spikes are visible in the anoma-
lous case, whereas profiles with more rounded maxima
occur for normal dispersion. A possible interpretation is
suggested by fig. 7, where we saw that anomalous dis-
persion tends to create regions of higher gradient.
The shocking region in fig. 10 has a similar outline
to that of fig. 9(b) except that it is slightly narrower,
as expected from the coherence length discussion above.
However, the shocking region in fig. 11 is much nar-
rower, especially given the change in the ∆n scale. This
is because the multi-stepped nature of the refractive in-
dex leads to a correspondingly large range of coherence
lengths, with shorter ones corresponding to those span-
ning several steps. Since the shocking region in the multi-
step case has reduced in size by a factor of two or three,
a reasonable inference might be that the dominant co-
herence length results from interaction over two or three
refractive index steps.
In realistic media, the refractive index will vary
smoothly with frequency, and the group velocity can dif-
ferent from the phase velocity. We can approximate this
situation most simply using a refractive index gradient
−0.1 −0.05 0 0.05 0.1
∆ n3,5
−0.1 −0.05 0 0.05 0.1
FIG. 10: Carrier shocking for CW and pulsed (τ = 0.93̇)
cases with a single refractive index step at 4ω1; showing the
correlation in the region where LC . LSPM between energy
conservation failure (logarithmic right hand scale, dots) and
the LDD detected shocking distance (left hand scale, solid
line). The dotted line shows the LDD results for a CW field
(τ = ∞).
−0.05 0 0.05
∆ n1,3,5...
−0.05 0 0.05
FIG. 11: Carrier shocking for CW and pulsed (τ = 0.93̇) cases
with multiple refractive index steps; showing the correlation
in the region where LC . LSPM between energy conserva-
tion failure (logarithmic right hand scale, dots) and the LDD
detected shocking distance (left hand scale, solid line). Note
the narrower range of ∆n as compared to the previous two
graphs. The dotted line shows the LDD results for a CW field
(τ = ∞).
rather than a series of steps. The results using the LDD
method in this case can be seen on fig. 12, where now
we also vary the strength of the nonlinearity. As in the
previous cases, we see a well defined shocking regime that
is asymmetric about the non-dispersive case.
Detailed examination shows that the curves in fig. 12
exhibit a marked similarity, and can be brought into
near perfect coincidence by applying the scaling χ(3) →
χ(3)/m, L → mL, and ∆n → ∆n/m. We also get a
comparable similarity for each of figs. 9,10,11, when sim-
ulations at nonlinear strengths of χ(3)E20 = 0.01 and
0.04 are added to the results. This demonstrates that
the character of the shocking is dominated by the sign of
inter-harmonic phase velocity differences, not by the local
group velocity dispersion at each harmonic. Thus anoma-
lous (normal) dispersion is primarily interesting because
7 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
it gives a negative (positive) refractive index shift be-
tween successive harmonics. We leave the complicated
(and far more subtle) effects of group velocity differences
or dispersion for later work.
FIG. 12: Shocking distance vs mismatch for weak refractive
index gradients δ, and a range of nonlinearities. This shows an
abrupt cut-off at positive ∆n, but a relatively gradual one at
negative ∆n. The initial pulses are identical to those in fig. 9.
No CW results are shown, since the dispersion experienced by
the field would be identical to the multi-stepped case shown
in fig. 11.
An important feature of all the results is the pro-
nounced asymmetry, with shocking persisting much fur-
ther into the anomalous dispersion regime; we have even
seen it for the case of a weakly parabolic refractive index.
We can quantify this asymmetry by considering lin-
ear and nonlinear contributions to the phase matching of
the harmonics (i.e. linear and nonlinear refractive index
shifts), as described at the beginning of section V. In the
simple single step case shown on fig. 9(b), the 1st–3rd
harmonic phase shift will dominate, because it applies to
the two most intense spectral components. We can calcu-
late the SPM-induced refractive index shift between the
fundamental and third harmonic with eqns. (18) and (19)
to be ∆nNL:1−∆nNL:3 = −0.011, since χ
(3)E2 = 0.02 at
the peak of the carrier oscillations, the refractive index is
n20 = 2, and 2n2E
1 = 2 × 3χ
(3)E2/8n0. This is roughly
comparable to the offset of the shocking region, which
is centred at about ∆n1,3 ≃ −0.014. We cannot ex-
pect perfect agreement, since the calculations ignore the
role of higher harmonic generation and depletion of the
fundamental. Whilst a simple calculation is reasonably
successful in this single-step case, it cannot be applied
for a realistic medium – or indeed to the situation shown
in fig. 12. There, the effects of dispersion, higher har-
monic generation, and nonlinear refractive index shifts
are inextricably intertwined.
To summarize, we have demonstrated that shocking
is strongly dependent on the interplay between LC and
LSPM , with shorter LC ’s (increasing ∆n’s) decreasing
the likelihood of shocking. We have also deduced the
reasons for the strong asymmetry of the shocking region.
VII. APPLICATIONS AND EXPERIMENTS
In sections I and IV, we discussed how imminent car-
rier shocking might be recognised computationally. In
considering whether shocking might be detectable exper-
imentally, we must now decide how it might manifest
itself in the laboratory. A mathematical discontinuity is
clearly not a physical possibility, even if Rosen [2] did
manage to accommodate it theoretically, albiet at the
expense of energy conservation. In practice, the increas-
ing field gradients (and spectral broadening) that precede
a shock will inevitably engender new physical processes
that will limit the steepening. Indeed, we have already
seen this happening in the previous two sections, where
dispersion has been seen to frustrate the self-steepening
process; the next barrier would be the time-scale of the
nonlinear response.
FIG. 13: Numerically predicted shocking distances Lshock
in fused silica. These are the MOC predicted distances for
the waveform, assuming the dispersion was (abruptly) ne-
glected, the shortest distance shown is about 1.9µm. The
damage threshold for fused silica was taken to be Pthreshold =
50TW/cm2.
Although there is no question of a mathematical dis-
continuity being observed in an experiment, the recent
advances in the measurment of optical pulse profiles (see
e.g. [13]), suggest that it might be possible to observe
carrier steepening. Indeed, it has already been pre-
dicted [4] that noticeable steepening effects could occur
for pulses propagating in fused silica. Unfortunately, in
our own simulations of this process, the LDD shock de-
tection was not triggered at any realistic pulse intensity.
Although visible steeping did occur, at no point was a
seriously distorted waveform approached, and the effects
would have been milder if we had included the finite re-
sponse time of the nonlinearity. The results presented
in fig. 13 show the steepest gradient recorded as a func-
tion of distance and pulse intensity for these simulations.
The third harmonic coherence length for these parame-
ters is about 7µm, and we can attribute the regular vari-
ation with distance seen in the figure to the third har-
monic component aligning with successive oscillations of
the fundamental as the waves move across each other.
8 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
The key reason why incipient shocking was not detected
is because the pulse frequency lies in a region of weak
normal dispersion, whereas we have shown in section V
that a region of weak anomalous dispersion (assuming it
could be achieved) would be favourable.
The current interest in media with tailored dispersion
characteristics [14, 15, 16] raises the interesting question
of whether it might be possible to engineer a material
that maximize self-steepening. A major stumbling block
would be the need for control over many harmonic orders;
however, for a narrow-band pulse, only the dispersion
characteristics close to the harmonics would be relevant,
which might perhaps make the technical challenge less
formidable.
The recognition that dispersion control can enhance
(or reduce) carrier self-steepening suggests other appli-
cations if we widen our horizons to encompass the more
general idea of carrier shaping. In this case, we would
exploit both nonlinearity and dispersion control to opti-
mize the shape of the carrier oscillations for a particular
experiment. Applications such as high harmonic genera-
tion (e.g. [17]) might well benefit from suitably designed
carrier wave modulation.
VIII. CONCLUSIONS
In this paper we have investigated carrier shock for-
mation, developed criteria for detecting its onset in nu-
merical simulations, and shown how it is influenced by
a range of parameters, particularly dispersion. We have
also obtained remarkable agreement between numerical
simulations and theoretical predictions of the shocking
distance in the dispersionless limit, and shown that the
process is sensitive to both CEP and pulse duration.
Although we have confirmed that shocking occurs in a
narrow parameter range, this is far from being the whole
picture. In particular, there is a distinct asymmetry be-
tween the anomalous and normal dispersion regimes. The
former leads to shocking signatures such as the appear-
ance of narrow spikes as the higher harmonics interfere
on the steepening part of the pulse profile. In contrast,
normal dispersion creates no such features, and the pulse
profiles have a rather blunt appearance. The asymmetry
arises from the effect of the nonlinear refractive index on
the dispersion induced phase mismatch.
The conclusion to be drawn from our results is clear: if
they could be engineered, materials with wide regions of
weakly anomalous dispersion are much better candidates
for generating steep, shock-like field profiles than those
(such as silica) with a weak normal dispersion. Detecting
incipient shock formation in materials like fused silica is
likely to be a near impossible task, given the constraints
imposed by their damage thresholds.
[1] F. DeMartini, C.H. Townes, T.K. Gustafson, P.L. Kelly,
Phys. Rev. 164, 312 (1967).
[2] G. Rosen, Phys. Rev. 139, A539 (1965).
[3] R.G. Flesch, A. Pushkarev, and J.V. Moloney, Phys. Rev.
Lett. 76, 2488 (1996).
[4] L. Gilles, J.V. Moloney, and L. Vázquez, Phys. Rev. E
60, 1051 (1999).
[5] S.Trillo, G. Millot, C. Conti, ECLEO 2006, paper EB1-
05-Mon.
[6] A. M. Kamchatnov, R. A. Kraenkel, B. A. Umarov,
Physi. Rev. E66, 036609 (2002).
[7] J.C.A. Tyrrell, P. Kinsler, G.H.C. New, J.Mod.Opt. 52,
973 (2005).
[8] B. Fornberg, “A Practical Guide to Pseudospectral Meth-
ods”, Cambridge University Press (1996).
[9] L. Gilles, S.C. Hagness, and L. Vázquez, J. Comp.Phys.
161 379 (2000).
[10] P. Kinsler, S.B.P. Radnor, G.H.C. New, Phys. Rev. A.
72, 063807 (2005).
[11] “Numerical Recipes in C”, W.H. Press et al., (CUP,
1992).
[12] H.-G.Pagendarm, B.Seitz, in “Scientific Visualization
- Advanced Software Techniques” (EllisHorwood 1993),
p159.
[13] E. Goulielmakis, M. Uiberacker, R. Kienberger, A. Bal-
tuska, V. Yakovlev, A. Scrinzi, Th. Westerwalbesloh, U.
Kleineberg, U. Heinzmann, M. Drescher, F. Krausz Sci-
ence 305, 1267 (2004).
[14] Photonic Bandgap Materials, (Ed: C. M. Soukoulis),
(Kluwer, Dordrecht, The Netherlands 1996).
[15] Photonic Crystals, (Eds: J. D. Joannopoulos, R. D.
Meade, J. N. Winn), (Princeton Academic Press, Prince-
ton, NJ 1995).
[16] J.K. Ranka, R.S. Windeler, A.J. Stentz, Opt. Lett. 25,
25-27 (2000).
[17] I.P. Christov, M.M. Murnane, H.C. Kapteyn, Phys. Rev.
Lett. 78, 1251 (1997).
9 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org

  Carrier wave shocking is studied using the Pseudo-Spectral Spatial Domain
(PSSD) technique. We describe the shock detection diagnostics necessary for
this numerical study, and verify them against theoretical shocking predictions
for the dispersionless case. These predictions show Carrier Envelope Phase
(CEP) and pulse bandwidth sensitivity in the single-cycle regime. The flexible
dispersion management offered by PSSD enables us to independently control the
linear and nonlinear dispersion. Customized dispersion profiles allow us to
analyze the development of both carrier self-steepening and shocks. The results
exhibit a marked asymmetry between normal and anomalous dispersion, both in the
limits of the shocking regime and in the (near) shocked pulse waveforms.
Combining these insights, we offer some suggestions on how carrier shocking (or
at least extreme self-steepening) might be realised experimentally.

arXiv:0704.1212v1  [physics.optics]  10 Apr 2007
CSHOCK-II Dr.Paul.Kinsler@physics.org
Optical carrier wave shocking: detection and dispersion
P. Kinsler, S.B.P. Radnor, J.C.A. Tyrrell, and G.H.C. New
Blackett Laboratory, Imperial College London, Prince Consort Road, London SW7 2BW, United Kingdom.
(Dated: November 20, 2018)
Carrier wave shocking is studied using the Pseudo-Spectral Spatial Domain (PSSD) technique.
We describe the shock detection diagnostics necessary for this numerical study, and verify them
against theoretical shocking predictions for the dispersionless case. These predictions show Carrier
Envelope Phase (CEP) and pulse bandwidth sensitivity in the single-cycle regime. The flexible
dispersion management offered by PSSD enables us to independently control the linear and nonlinear
dispersion. Customized dispersion profiles allow us to analyze the development of both carrier self-
steepening and shocks. The results exhibit a marked asymmetry between normal and anomalous
dispersion, both in the limits of the shocking regime and in the (near) shocked pulse waveforms.
Combining these insights, we offer some suggestions on how carrier shocking (or at least extreme
self-steepening) might be realised experimentally.
I. INTRODUCTION
The self-steepening of an optical pulse envelope was
first studied by DeMartini et al. in 1967 [1], and is a
well-known phenomenon associated with self-phase mod-
ulation (SPM). Surprisingly however, the possibility of
self-steepening of the optical carrier wave1 was considered
even earlier in a 1965 paper by Rosen [2], who showed
that, for a third order χ(3) nonlinearity (and under suit-
able conditions), a shock (or field discontinuity) develops
in a finite distance. This latter phenomenon received lit-
tle attention for more than 30 years, until it was revisited
in the 1990s by Moloney et al. [3, 4], who performed Fi-
nite Difference Time Domain (FDTD) simulations of the
process. Three dimensional FDTD simulations of carrier
shocking have recently been performed by Trillo et al.
In the present paper, we investigate carrier wave shock-
ing in χ(3) nonlinear materials where SPM is accompa-
nied by the generation of (odd) higher harmonics. For
the dispersionless case we have generalized earlier pre-
dictions based on the method of characteristics (MOC)
to allow for arbitrary initial waveforms. This allows us,
for example, to predict the carrier envelope phase (CEP)
sensitivity of the shocking distance for optical pulses, as
well as the dependence on pulse length.
Our primary interest in this paper is the effect of linear
dispersion on carrier shock formation. Since no analytic
solutions exist in this case, we are forced to rely on nu-
merical simulations. In any discussion of carrier wave
shocking, it is important to distinguish between three
related concepts: the physical system, the mathemati-
cal model, and the numerical model. The mathematical
model is an approximation to the physical system, while
∗Electronic address: Dr.Paul.Kinsler@physics.org
1 Note that in this paper we use the term “carrier” to denote
all the oscillations of the field, and do not use its other sense,
i.e. that of fixed-frequency oscillations as used in a envelope and
carrier representation of a pulse.
the numerical model is an approximation to the mathe-
matical model. Actual discontinuities occur only in the
mathematical model, and it is these to which the idea of
shock formation refers. In the physical system, disconti-
nuities are prevented by phenomena not included in the
mathematics, while numerical codes inevitably fail as a
mathematical discontinuity is approached. The indicator
of an imminent “shock” is the rapidly increasing gradi-
ent (steepening) of the optical carrier. In the numerical
model, we see numerical symptoms generated by extreme
self-steepening, and these correlate with the onset of the
discontinuities in the mathematical model. Under these
circumstances, it has been necessary for us to develop
a quantitative numerical test of “shock formation”. We
have found the most satisfactory diagnostic to be “Local
Discontinuity Detection” (LDD), which gives results that
are in good agreement with theoretical MOC predictions
based on the mathematical model. LDD provides a clear
numerical measure of the rapid steepening that precedes
the appearance of a discontinuity in the mathematical
representation. We continue to use the term “shock” to
describe the situation where a mathematical shock is im-
minent.
In our simulations, we exploit the flexibility in disper-
sion management offered by the Pseudo-Spectral Spatial
Domain (PSSD) technique [7] to study carrier shock for-
mation for a range of simple dispersion profiles, and de-
termine the degree of phase mismatch that carrier shock-
ing can tolerate. As expected, we find that shocking oc-
curs when the nonlinearity dominates the linear disper-
sion. However, it emerges that the process is asymmet-
ric, with anomalous dispersion being far more conducive
to shocking than normal dispersion. Hence, the fact that
the LDD scheme does not detect carrier shocks in simula-
tions involving (normally dispersive) fused silica, even at
powers equal to its damage threshold, is neither surpris-
ing nor necessarily discouraging. After all, anomalously
dispersive materials could potentially be engineered. Fur-
ther, our results also relate to how one might perform
carrier shaping (as opposed to carrier steepening), a pro-
cess that has some interesting applications.
1 Kinsler-RTN-2007
http://arxiv.org/abs/0704.1212v1
mailto:Dr.Paul.Kinsler@physics.org
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
After briefly describing our simulation methods in
section II, we consider dispersionless shocking and the
method of characteristics in III, and our LDD shock de-
tection scheme in IV. Then we discuss the effect of disper-
sion on carrier shocking in V, followed by our numerical
results in VI. In VII we consider the potential relevance
to experimental detection of carrier steepening and/or
shocks. Finally, in VIII, we present our conclusions.
II. SIMULATION METHODS
The PSSD method [7, 8] offers significant advantages
over the traditional FDTD and Pseudospectral Time-
Domain (PSTD) [9] techniques for modeling the propa-
gation and interaction of few-cycle pulses. Run times are
generally faster, and PSSD also offers far greater flexi-
bility in the handling of dispersion. Whereas FDTD and
PSTD [9] propagate fields E(z), H(z) forward in time,
PSSD propagates fields E(t), H(t) forward in space. It
is important to keep this difference in mind when com-
paring our results to those in [9]. Under PSSD, the en-
tire time-history (and therefore frequency content) of the
pulse is known at any point in space, so arbitrary disper-
sion incurs no extra computational penalty. In contrast,
the FDTD or PSTD approaches use convolutions incor-
porating time-response models for dispersion.
We apply the PSSD algorithm to two representations
of the field and source-free Maxwell’s equations in non-
magnetic media; the first uses the E and H fields, and
the second the directional fields ~G±(t) = αr(t)∗ ~Ex(t)±
βr ~Hy(t) [10]. Here the αr, βr include the (linear) per-
mittivity and permeability of the material (i.e. ǫ(t), µ).
These G± fields enable us to rewrite Maxwell’s equations,
and efficiently separate out the relevant forward-going
part of the field.
For an instantaneous χ(3) nonlinearity, the equations
for E and H in the 1D (plane wave) limit are
dHy(t; z)
ǫ0ǫr(t) ∗ Ex(t; z) + χ
(3)Ex(t; z)
dEx(t; z)
[µ0Hy(t; z)] . (2)
The G± field simulations usually assume G− = 0, and as
a result contain only forward traveling components. The
forward-only wave equation for G+ is
dG+(t; z)
βrαr(t) ∗G
+(t; z) + βrχ
(3)E(t; z)3
where it is most straightforward to calculate the non-
linear term by reconstructing E(t; z) from G+(t; z) in
the frequency domain using Ẽ(ω; z) = G̃+(ω; z)/2α̃r(ω),
since G− = 0. Notice the similarity between eqns. (1)
and (3), but that eqn. (3) propagates the field in a single
first order equation, rather than two.
Typical array sizes used in pulse simulations were
N = 214 covering a time window T = 200fs, and (spatial)
propagation steps were dz = 0.4cT/N ≈ 0.9nm. We en-
sured the stability of our integration using Orszag’s 2/3
rule [11], which involves setting the upper part of the
spectral range to zero. It is worth noting that changing
this cut-off, either by adjusting its position, or by using
a smoothed (rather than step-like) filter, made little dif-
ference to test simulations. The pulse profile used as an
initial condition was
E(t) = E0 sin(ω1t+ φ) sech(0.28ω1t/τ), (4)
where our standard parameters were ω1 = 2.356 ×
1015rad/s (i.e. λ = 800nm) with τ = 0.93̇. Such pulses
are rather short (since the number of cycles inside the
intensity FWHM is τ), but in fact the shocking distance
is only weakly dependent on the pulse width, with signif-
icant variation only appearing for pulses of a few cycles
or less.
We also performed CW simulations, for which we mod-
eled just a single cycle of the carrier. assisted by the pe-
riodic nature of the discrete Fourier transform. For these
we used array sizes of N = 210, and the time window was
set by the period of the field oscillations.
Our default value of nonlinear strength was χ(3)E20 =
0.02, which is comparable to that in fused silica at an in-
tensity of 0.7×1014 W/cm2; our χ(3) parameter is equiv-
alent to η in Rosen [2], and a in Gilles et al.[4]. We use an
instantaneous nonlinearity, since our primary interest is
linear dispersion, and that is a far more significant effect
than the nonlinear time response.
III. CARRIER WAVE SHOCKING
As an introduction to the process of carrier wave shock-
ing, fig. 1 shows the profile of a pulse propagating in a
dispersionless χ(3) medium, just before a shock occurs.
The nonlinearity gives rise to a nonlinear index of refrac-
tion n2E
2, the effect of which is to increase the effective
refractive index in the more intense regions of the pro-
file. This reduces the phase velocity at the peak of each
oscillation with respect to the rest of the waveform, and
causes the slope on the trailing edges to increase dramat-
ically.
The effects seen in fig. 1 are associated with the gen-
eration of third and higher harmonics, although the har-
monic components are not particularly strong even when
a shock is about to occur. As an example, fig. 2 shows
how the harmonics build up as shocking is approached.
The ω4 scaling of the intensity spectrum in the figure ex-
aggerates the contribution of the higher orders, and has
been chosen for illustrative purposes. Notice that the
profile becomes nearly flat (i.e. the spectrum falls off as
the 4th power of the frequency) just before shocking is
registered at around 4.3µm.
If we write
E = A(t) [sin(ω1t) + γ cos(3ω1t+ ψ)] , (5)
2 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
5 10 15 20
time(fs)
FIG. 1: The profile of a few-cycle optical pulse just prior to
shocking in the dispersionless limit. The larger oscillations
in the centre of the pulse undergo more self-steepening than
those in the wings. The standard pulse parameters were used.
FIG. 2: Development of the heights of the scaled harmonic
peaks as a pulse approaches the (LDD) shocking distance of
4.3µm in the dispersionless case. Each line perpendicular to
the ω/ω1 axis corresponds to the (scaled) contribution from
that spectral peak. Note that the viewpoint has been rotated
so that the contribution from the fundamental is to the right.
The initial pulse contained about 33 cycles (τ = 33).
we find that an appropriate choice of A(t), with γ =
0.1 and ψ = 0, gives us a passable match to fig. 1.
Note that this choice of ψ corresponds to the phase of
third harmonic generation (THG) under index matched
conditions, i.e. where n0 = n(ω1) = n(ω3).
Rosen’s original paper [2] used the MOC to predict
the formation of a value discontinuity in the field at cer-
tain points within the profile. If the displacement of the
dispersionless medium is written
D = ǫ0
E + χ(1)E + χ(3)E3
, (6)
he showed that the wave equation for E is
1 + χ(1)
) ∂2E
+ χ(3)
, (7)
and that the associated equation governing the charac-
teristic lines of E is
+ v(E)
= 0. (8)
Here, the velocity v(E) is given by
v(E) =
(ǫr + 3χ(3)E2)1/2
, (9)
where ǫr = 1 + χ
(1) = n0
2 is the (relative) dielectric
constant and n0 the linear refractive index.
v − dv
FIG. 3: Method of Characteristics. Two points A and B
on the field profile, separated initially by a time difference dt
travel at different speeds v − dv and v, and meet at point C.
Using eqn. (9) along with the construction shown in
fig. 3, we can derive a simple formula for the distance to
shocking. The figure shows two characteristics AC and
BC, originating from points A and B, and converging
towards a shock at C after a distance of L. The intensity
at A is higher than at B, so the speed associated with
AC (represented by its gradient) is lower than that of
BC. This means that at C, the field has two values, and
a discontinuity has formed. From the geometry of the
figure, it is easy to show that
(c/n0)
, (10)
where t, v = c/n0, and L = vt are respectively time,
speed and distance.
On the other hand, differentiating eqn. (9) leads to
3cχ(3)
n2o + 3χ
(3)E2
, (11)
and this combined with eqn. (10) yields
1 + 3χ(3)E2/n2o
3χ(3)(−dE2/dt)
1 + 8n2E2/no
(−dE2/dt)
, (13)
where n2 = 3χ
(3)/8no is the material parameter deter-
mining the intensity induced refractive index shift as
2. For a given profile, a shock will occur first at
the point where −dE2/dt reaches its negative extremum.
3 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
We can therefore define the shocking distance as
Lshock =
1 + 8n2E2/no
(−dE2/dt)
(−dE2/dt)
for 8n2E
2/no ≪ 1.(15)
This formula is more general than that of either Rosen
[2] or Gilles et al. [4]. Notice in particular that the
parameter that controls the shock behaviour is not the
gradient of the field (dE/dt), but that of the field squared
(dE2/dt).
For a sinusoidal initial waveform E = E0 sin(ω1t+ φ),
it is easy to show that shocks form at
ω1t = −π/4 + jπ − φ, (16)
where j is an integer, and that the shocking distance is
Lshock =
3ω1χ(3)E
4ω1n2E
. (17)
The approach is readily extended to pulsed waveforms
of the type defined in eqn. (4). Analytical results can
be derived in the new situation in which eqn. (16) be-
comes a transcendental equation. However, the results
are cumbersome and it is simpler to scan the profile nu-
merically to determine the shocking parameters. It turns
out that the shocking distance for short pulses exhibits
an interesting sensitivity to the carrier envelope phase φ.
Whilst in the CW case all locations defined by eqn. (16)
were equivalent, for pulses, the one nearest the peak of
the envelope has a shorter Lshock than the others.
A set of results is displayed in fig. 4 where the shock-
ing distance is plotted as a function of the carrier phase
φ for sech profiles with different pulse widths τ . The
dotted line is for the case of a very broad pulse for which
Lshock is given by eqn. (17). The sharp peaks mark a
curve crossings where the shock location switches from
one point on the sinusoid to another. Notice that the
range of shocking distances increases for shorter pulses,
and that, unsurprisingly, the lowest values occur when φ
is around π/2 i.e. when the pulse has a cosine form.
One very important point to note from fig. 4 is that for
pulses containing more than a few cycles, the dependence
of the shocking distance on pulse width is very weak. We
will see this message repeated later in section V, with
shock regions being similar for both single cycle (τ = 1)
pulses and CW (τ = ∞) fields.
IV. SHOCK DETECTION
Since optical shock formation is directly associated
with regions of increasingly steep field gradient, any nu-
merical scheme, however sophisticated, is bound to fail at
some point in the process. We therefore want to recognize
when a shock is imminent, not only to avoid numerical
FIG. 4: MOC shocking distances as a function of phase,
for pulses as in eqn. (4), allowing for different pulse lengths.
We include the 1/2 cycle τ = 0.56 results to emphasize the
trend; but even for a long pulse (τ = 28), a small peak can
be seen just below φ/π = 0.01. The peak position is weakly
τ dependent. The + signs denote LDD shocking distances
obtained from simulations of the τ = 0.93̇ case.
problems, but as a means to estimate the distance at
which a discontinuity would occur in the mathematical
model.
One obvious symptom of impending numerical failure
is loss of energy conservation [7]. However, this does not
give an indication at the first instance of a shock form-
ing, but rather signals the accumulated effect of multiple
small numerical failures from many shocked regions.
A more physical strategy is to search for regions where
the field gradient dE/dt is large and increasing rapidly,
and to use this to predict the shocking point. A useful
variant, suggested by the MOC calculation in the previ-
ous section, is to use the value of −dE2/dt instead of the
gradient.
Overall, we find that the best method is Local Discon-
tinuity Detection (LDD), which is similar to techniques
used in other fields (see e.g. [12]). As the shock regime is
approached, narrow shoulders with associated points of
inflection appear within the regions of rapidly increasing
gradient. The procedure is therefore to scan the field pro-
file for the maximum gradient (of either E or E2) and, if
it occurs near a point of inflection, an incipient shock is
registered. An example of a pulse that has just triggered
the LDD diagnostic is shown in fig 5.
The LDD method requires two parameters. The first
determines the time scale used in determining whether a
point of inflection exists. For this, we pick the scale set
by our temporal grid and insist that the field gradients
calculated at three adjacent grid points have opposing
signs: either up-down-up, or down-up-down. The sec-
ond determines the maximum range allowed between the
maximum gradient and the point of inflection, and our
default value for this was 10 grid points. In our simula-
tions, we see that the position of the first detected shock
depends only weakly on this range. We can easily mini-
mize the small sensitivity to these parameters by holding
4 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
them fixed throughout any given set of simulations. As a
result, we have found LDD to be a sensitive and reliable
method of shock detection.
FIG. 5: The profile of a few-cycle optical pulse at the point
of shocking. The larger oscillations in the centre of the pulse
undergo more self-steepening than those in the wings. (a) The
whole pulse, with the LDD carrier shock circled (◦). (b) An
expanded view of the shock region, with its point of inflection
(at ∼ 9.3fs) being very close to the steepest gradient (at ∼
9.1fs). The standard pulse parameters were used.
Eqn. (17) predicts that the shocking distance should
increase linearly with the refractive index n0. We use
this to test the LDD diagnostic in fig. 6, where the ana-
lytical formula is compared with the results of numerical
simulations. This figure shows close agreement between
prediction (using eqn. (15)) and simulation, where the
approximation causes the MOC prediction to be reduced
by less than 0.1µm at n0 = 1. We see similar agreement
between simulations using the LDD method and the car-
rier phase sensitivity shown in fig. 4. The presence of
small systematic differences (as on e.g. fig. 4, 6) can be
easily understood, since the LDD diagnostic is (strictly
speaking) a test of the numerics, and is not a direct test
for the presence of a physical shock or mathematical dis-
continuity.
V. THE EFFECT OF DISPERSION ON
SHOCKING
The primary purpose of this paper is to understand
the principles of carrier shock formation in the presence
of dispersion. Some simple ideas about the role of dis-
persion can be understood from eqn.(5) using the insight
FIG. 6: The LDD shocking distance as a function of refractive
index, comparing the approximate MOC prediction from eqn.
(15) to PSSD simulations (+) for pulses with τ = 0.93̇. Other
PSSD simulation results from independent codes give very
similar results to those shown on the graph.
from eqns.(10) – (15), i.e. that the rate of change of E2
is the critical factor in shock development, rather than
that of E. If only the fundamental and third harmonic
components are considered (as in eqn.(5)), the effective
refractive indices are
∆nNL:1 = n1 + n2
1 + 3γ2
I1, (18)
∆nNL:3 = n3 + n2
3 + γ2
I1. (19)
Evidently, the relative phase velocity of the two waves is
affected by both linear and nonlinear dispersion, so the
phase ψ in eqn.(5) will vary accordingly as the pulse prop-
agates. Fig. 7, which shows how dE2/dt varies with ψ,
suggests that shocking is likely to be exacerbated when ψ
is small and positive, but moderated when ψ is negative.
Broadly speaking, the former case will be promoted by
anomalous dispersion and the latter by normal disper-
sion; however, the process is clearly complicated, since it
involves time-dependent phase shifts between the waves
and the interplay of linear and nonlinear dispersion. Of
course, since carrier shocking relies on the establishment
and maintenance of specific phase relationships between
a set of harmonics, it must be expected that strong dis-
persion of either sign will disrupt the shock formation
process. On the other hand, the simple argument that
has been offered suggests that shocking may be tolerant
to a degree of anomalous dispersion, but not to a simi-
lar amount of normal dispersion. In general, therefore, a
graph of shocking signature versus refractive index mis-
match might be expected to exhibit a shock region where
nonlinearity dominates dispersion (displaced in the direc-
tion of anomalous dispersion), surrounded by a shock-
free region where dispersion dominates the nonlinearity.
Moreover, if a dominant coherence length LC can be de-
fined, it is reasonable to expect shocking to occur when
this exceeds the characteristic SPM length (LSPM ). As
we shall see, all these features are borne out by the nu-
merical results.
5 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
FIG. 7: The effect of different time-lags between fundamental
and third harmonic on the steepness of the pulse. The field
profile is as in eqn. (5), with A(t) = A0 and γ = 1/16. (a)
Scaled maximum values of Mj = Max(dE
j/dt) as a function
of third harmonic phase offset ψ. We can see that the maxima
of dE/dt and dE2/dt occur at different ψ. (b) Position on the
pulse θj = ω1tj of the maxima plotted in the top frame. Note
that the kinks in (a) and (b) occur for offsets that give the
longest shocking distances, not the shortest.
In our numerical simulations, we make extensive use
of model dispersion curves, which enable refractive in-
dex differences between harmonics to be freely controlled,
and lead to uncomplicated boundaries between the shock
and no-shock regimes. The dispersion profiles shown
in fig. 8 duly contain either refractive index steps ∆n
at the midpoints between successive odd harmonics, or
smooth gradients δ that provide a similar net change.
In the simplest option, there is a single step at 2ω1,
in which case the dominant coherence length is clearly
LC = π/ |k3 − 3k1| = πc/3ω1 |n3 − n1|; we have also
tried a single step at 4ω1, which has a shorter LC . In
all cases, when n increases with frequency, the situation
corresponds broadly to normal dispersion, while decreas-
ing n corresponds to anomalous dispersion. The step
size (or gradient) is chosen to be comparable to values in
fused silica at ω1 = 1.5rad/fs, where the refractive index
differences between the lower harmonics are ∆n1,3 ≈ 0.06
and ∆n3,5 ≈ 0.12[4]. While we consider the case of fused
silica itself in Section VII, the results based on the model
dispersion characteristics are invaluable for understand-
2   n∆
ω1 ω3 ω5 ω7
FIG. 8: Types of refractive index profile used. The solid lines
show a single refractive index step midway between ω1 and
ω3, so that the fundamental is always phase mismatched from
its higher harmonics, although they remain perfectly phase
matched with each other. The dashed lines show multiple re-
fractive index steps, each of the same size, and always midway
between subsequent harmonics. The dotted lines show a lin-
ear refractive index gradient which gives the same mismatch
between subsequent harmonics as the multi-stepped case. We
do not show lines for the case of a single step at 4ω1 to avoid
cluttering the figure.
ing the essential principles of carrier shock formation in
dispersive media.
VI. RESULTS AND DISCUSSION
We will now analyze our numerical simulations of car-
rier shocking in the presence of dispersion on the basis of
the principles discussed above. Results are included for
both a CW wave and for a single-cycle pulse. The CW
case gives slightly wider shocking regions, but the differ-
ences are minor. This is because the effect of the pulse
envelope on the field amplitudes and gradients of the cen-
tral carrier oscillation is small, except when considering
sub-cycle pulses. In the results we present, the energy
conservation and LDD measures reveal the presence of
sharp boundaries between the shocking and non-shocking
regions.
Results for single refractive index steps at 2ω1 and 4ω1
are presented in figs. 9 and 10 respectively, whereas fig.
11 has a step midway between all harmonics. In both
single stepped cases, a useful coherence length can be
defined on the basis of the mismatch between the fun-
damental and the harmonic just above the step. In the
multi-stepped case, there is no easy way to define a dom-
inant coherence length.
Fig. 9(a) shows the pulse profiles in the anomalous
(negative) and normal (positive) single step cases, for the
smallest step size at which shocking did not occur. The
first obvious difference between the profiles is the oppo-
site walk-off direction of the third harmonic, although the
third harmonic contribution is hard to see in the anoma-
lous case. The second is that the pulse profiles exhibit
6 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
−20 −15 −10 −5 0
time(fs)
−0.1 −0.05 0 0.05 0.1
∆ n1,3
−0.1 −0.05 0 0.05 0.1
(b) X
FIG. 9: Carrier shocking for CW and pulsed (τ = 0.93̇) cases
with a single refractive index step at 2ω1. Frame (a) compares
chirped pulses at a 40µm propagation distance from the non-
shock region immediately outside of the shocked region; the
upper curve (X) is for the negative step, the lower (Y) for
the positive step. Frame (b) shows the correlation in the
region where LC . LSPM between energy conservation failure
(logarithmic right hand scale, dots) and the LDD detected
shocking distance (left hand scale, solid line). The dotted
line shows the LDD results for a CW field (τ = ∞).
distinctly different characteristics according to the sign of
the dispersion. Narrow spikes are visible in the anoma-
lous case, whereas profiles with more rounded maxima
occur for normal dispersion. A possible interpretation is
suggested by fig. 7, where we saw that anomalous dis-
persion tends to create regions of higher gradient.
The shocking region in fig. 10 has a similar outline
to that of fig. 9(b) except that it is slightly narrower,
as expected from the coherence length discussion above.
However, the shocking region in fig. 11 is much nar-
rower, especially given the change in the ∆n scale. This
is because the multi-stepped nature of the refractive in-
dex leads to a correspondingly large range of coherence
lengths, with shorter ones corresponding to those span-
ning several steps. Since the shocking region in the multi-
step case has reduced in size by a factor of two or three,
a reasonable inference might be that the dominant co-
herence length results from interaction over two or three
refractive index steps.
In realistic media, the refractive index will vary
smoothly with frequency, and the group velocity can dif-
ferent from the phase velocity. We can approximate this
situation most simply using a refractive index gradient
−0.1 −0.05 0 0.05 0.1
∆ n3,5
−0.1 −0.05 0 0.05 0.1
FIG. 10: Carrier shocking for CW and pulsed (τ = 0.93̇)
cases with a single refractive index step at 4ω1; showing the
correlation in the region where LC . LSPM between energy
conservation failure (logarithmic right hand scale, dots) and
the LDD detected shocking distance (left hand scale, solid
line). The dotted line shows the LDD results for a CW field
(τ = ∞).
−0.05 0 0.05
∆ n1,3,5...
−0.05 0 0.05
FIG. 11: Carrier shocking for CW and pulsed (τ = 0.93̇) cases
with multiple refractive index steps; showing the correlation
in the region where LC . LSPM between energy conserva-
tion failure (logarithmic right hand scale, dots) and the LDD
detected shocking distance (left hand scale, solid line). Note
the narrower range of ∆n as compared to the previous two
graphs. The dotted line shows the LDD results for a CW field
(τ = ∞).
rather than a series of steps. The results using the LDD
method in this case can be seen on fig. 12, where now
we also vary the strength of the nonlinearity. As in the
previous cases, we see a well defined shocking regime that
is asymmetric about the non-dispersive case.
Detailed examination shows that the curves in fig. 12
exhibit a marked similarity, and can be brought into
near perfect coincidence by applying the scaling χ(3) →
χ(3)/m, L → mL, and ∆n → ∆n/m. We also get a
comparable similarity for each of figs. 9,10,11, when sim-
ulations at nonlinear strengths of χ(3)E20 = 0.01 and
0.04 are added to the results. This demonstrates that
the character of the shocking is dominated by the sign of
inter-harmonic phase velocity differences, not by the local
group velocity dispersion at each harmonic. Thus anoma-
lous (normal) dispersion is primarily interesting because
7 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
it gives a negative (positive) refractive index shift be-
tween successive harmonics. We leave the complicated
(and far more subtle) effects of group velocity differences
or dispersion for later work.
FIG. 12: Shocking distance vs mismatch for weak refractive
index gradients δ, and a range of nonlinearities. This shows an
abrupt cut-off at positive ∆n, but a relatively gradual one at
negative ∆n. The initial pulses are identical to those in fig. 9.
No CW results are shown, since the dispersion experienced by
the field would be identical to the multi-stepped case shown
in fig. 11.
An important feature of all the results is the pro-
nounced asymmetry, with shocking persisting much fur-
ther into the anomalous dispersion regime; we have even
seen it for the case of a weakly parabolic refractive index.
We can quantify this asymmetry by considering lin-
ear and nonlinear contributions to the phase matching of
the harmonics (i.e. linear and nonlinear refractive index
shifts), as described at the beginning of section V. In the
simple single step case shown on fig. 9(b), the 1st–3rd
harmonic phase shift will dominate, because it applies to
the two most intense spectral components. We can calcu-
late the SPM-induced refractive index shift between the
fundamental and third harmonic with eqns. (18) and (19)
to be ∆nNL:1−∆nNL:3 = −0.011, since χ
(3)E2 = 0.02 at
the peak of the carrier oscillations, the refractive index is
n20 = 2, and 2n2E
1 = 2 × 3χ
(3)E2/8n0. This is roughly
comparable to the offset of the shocking region, which
is centred at about ∆n1,3 ≃ −0.014. We cannot ex-
pect perfect agreement, since the calculations ignore the
role of higher harmonic generation and depletion of the
fundamental. Whilst a simple calculation is reasonably
successful in this single-step case, it cannot be applied
for a realistic medium – or indeed to the situation shown
in fig. 12. There, the effects of dispersion, higher har-
monic generation, and nonlinear refractive index shifts
are inextricably intertwined.
To summarize, we have demonstrated that shocking
is strongly dependent on the interplay between LC and
LSPM , with shorter LC ’s (increasing ∆n’s) decreasing
the likelihood of shocking. We have also deduced the
reasons for the strong asymmetry of the shocking region.
VII. APPLICATIONS AND EXPERIMENTS
In sections I and IV, we discussed how imminent car-
rier shocking might be recognised computationally. In
considering whether shocking might be detectable exper-
imentally, we must now decide how it might manifest
itself in the laboratory. A mathematical discontinuity is
clearly not a physical possibility, even if Rosen [2] did
manage to accommodate it theoretically, albiet at the
expense of energy conservation. In practice, the increas-
ing field gradients (and spectral broadening) that precede
a shock will inevitably engender new physical processes
that will limit the steepening. Indeed, we have already
seen this happening in the previous two sections, where
dispersion has been seen to frustrate the self-steepening
process; the next barrier would be the time-scale of the
nonlinear response.
FIG. 13: Numerically predicted shocking distances Lshock
in fused silica. These are the MOC predicted distances for
the waveform, assuming the dispersion was (abruptly) ne-
glected, the shortest distance shown is about 1.9µm. The
damage threshold for fused silica was taken to be Pthreshold =
50TW/cm2.
Although there is no question of a mathematical dis-
continuity being observed in an experiment, the recent
advances in the measurment of optical pulse profiles (see
e.g. [13]), suggest that it might be possible to observe
carrier steepening. Indeed, it has already been pre-
dicted [4] that noticeable steepening effects could occur
for pulses propagating in fused silica. Unfortunately, in
our own simulations of this process, the LDD shock de-
tection was not triggered at any realistic pulse intensity.
Although visible steeping did occur, at no point was a
seriously distorted waveform approached, and the effects
would have been milder if we had included the finite re-
sponse time of the nonlinearity. The results presented
in fig. 13 show the steepest gradient recorded as a func-
tion of distance and pulse intensity for these simulations.
The third harmonic coherence length for these parame-
ters is about 7µm, and we can attribute the regular vari-
ation with distance seen in the figure to the third har-
monic component aligning with successive oscillations of
the fundamental as the waves move across each other.
8 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
CSHOCK-II Dr.Paul.Kinsler@physics.org
The key reason why incipient shocking was not detected
is because the pulse frequency lies in a region of weak
normal dispersion, whereas we have shown in section V
that a region of weak anomalous dispersion (assuming it
could be achieved) would be favourable.
The current interest in media with tailored dispersion
characteristics [14, 15, 16] raises the interesting question
of whether it might be possible to engineer a material
that maximize self-steepening. A major stumbling block
would be the need for control over many harmonic orders;
however, for a narrow-band pulse, only the dispersion
characteristics close to the harmonics would be relevant,
which might perhaps make the technical challenge less
formidable.
The recognition that dispersion control can enhance
(or reduce) carrier self-steepening suggests other appli-
cations if we widen our horizons to encompass the more
general idea of carrier shaping. In this case, we would
exploit both nonlinearity and dispersion control to opti-
mize the shape of the carrier oscillations for a particular
experiment. Applications such as high harmonic genera-
tion (e.g. [17]) might well benefit from suitably designed
carrier wave modulation.
VIII. CONCLUSIONS
In this paper we have investigated carrier shock for-
mation, developed criteria for detecting its onset in nu-
merical simulations, and shown how it is influenced by
a range of parameters, particularly dispersion. We have
also obtained remarkable agreement between numerical
simulations and theoretical predictions of the shocking
distance in the dispersionless limit, and shown that the
process is sensitive to both CEP and pulse duration.
Although we have confirmed that shocking occurs in a
narrow parameter range, this is far from being the whole
picture. In particular, there is a distinct asymmetry be-
tween the anomalous and normal dispersion regimes. The
former leads to shocking signatures such as the appear-
ance of narrow spikes as the higher harmonics interfere
on the steepening part of the pulse profile. In contrast,
normal dispersion creates no such features, and the pulse
profiles have a rather blunt appearance. The asymmetry
arises from the effect of the nonlinear refractive index on
the dispersion induced phase mismatch.
The conclusion to be drawn from our results is clear: if
they could be engineered, materials with wide regions of
weakly anomalous dispersion are much better candidates
for generating steep, shock-like field profiles than those
(such as silica) with a weak normal dispersion. Detecting
incipient shock formation in materials like fused silica is
likely to be a near impossible task, given the constraints
imposed by their damage thresholds.
[1] F. DeMartini, C.H. Townes, T.K. Gustafson, P.L. Kelly,
Phys. Rev. 164, 312 (1967).
[2] G. Rosen, Phys. Rev. 139, A539 (1965).
[3] R.G. Flesch, A. Pushkarev, and J.V. Moloney, Phys. Rev.
Lett. 76, 2488 (1996).
[4] L. Gilles, J.V. Moloney, and L. Vázquez, Phys. Rev. E
60, 1051 (1999).
[5] S.Trillo, G. Millot, C. Conti, ECLEO 2006, paper EB1-
05-Mon.
[6] A. M. Kamchatnov, R. A. Kraenkel, B. A. Umarov,
Physi. Rev. E66, 036609 (2002).
[7] J.C.A. Tyrrell, P. Kinsler, G.H.C. New, J.Mod.Opt. 52,
973 (2005).
[8] B. Fornberg, “A Practical Guide to Pseudospectral Meth-
ods”, Cambridge University Press (1996).
[9] L. Gilles, S.C. Hagness, and L. Vázquez, J. Comp.Phys.
161 379 (2000).
[10] P. Kinsler, S.B.P. Radnor, G.H.C. New, Phys. Rev. A.
72, 063807 (2005).
[11] “Numerical Recipes in C”, W.H. Press et al., (CUP,
1992).
[12] H.-G.Pagendarm, B.Seitz, in “Scientific Visualization
- Advanced Software Techniques” (EllisHorwood 1993),
p159.
[13] E. Goulielmakis, M. Uiberacker, R. Kienberger, A. Bal-
tuska, V. Yakovlev, A. Scrinzi, Th. Westerwalbesloh, U.
Kleineberg, U. Heinzmann, M. Drescher, F. Krausz Sci-
ence 305, 1267 (2004).
[14] Photonic Bandgap Materials, (Ed: C. M. Soukoulis),
(Kluwer, Dordrecht, The Netherlands 1996).
[15] Photonic Crystals, (Eds: J. D. Joannopoulos, R. D.
Meade, J. N. Winn), (Princeton Academic Press, Prince-
ton, NJ 1995).
[16] J.K. Ranka, R.S. Windeler, A.J. Stentz, Opt. Lett. 25,
25-27 (2000).
[17] I.P. Christov, M.M. Murnane, H.C. Kapteyn, Phys. Rev.
Lett. 78, 1251 (1997).
9 Kinsler-RTN-2007
mailto:Dr.Paul.Kinsler@physics.org
