The Whitham approach to dispersive shocks in systems with cubic–quintic nonlinearities

By employing a rigorous approach based on the Whitham modulation theory, we investigate dispersive shock waves arising in a high-order nonlinear Schrödinger equation with competing cubic and quintic nonlinear responses. This model finds important applications in both nonlinear optics and Bose–Einstein condensates. Our theory predicts the formation of dispersive shocks with totally controllable properties, encompassing both steering and compression effects. Numerical simulations confirm these results perfectly. Quite remarkably, shock tuning can be achieved in the regime of a very small high order, i.e. quintic, nonlinearity.


Introduction
Dispersive shock waves (DSWs) are nonlinear phenomena characterized by the tendency of a wavepacket to form multivalued wavefronts, which are regularized by the appearance of nonstationary oscillating wave trains [1,2]. Multivalued regions are predicted by the hydrodynamic (dispersionless) approximation of the original model [3,4], which turns out to be a good description of the initial stage of propagation where dispersive effects are comparatively weak. The multivalued solutions in the hydrodynamic approximation come after the first (critical) time where a smooth initial profile develops an infinite steepness at some location, which defines the wave breaking (gradient catastrophe) point. After such wave breaking, the hydrodynamic approximation no longer describes the dynamics and an appropriate regularization incorporating dispersive effects must be used. In particular, the presence of dispersion, which becomes important close to the points of infinite steepness, is responsible for the spontaneous generation of expanding modulated periodic waves, which can be regarded as the natural oscillatory modes of nonlinear dispersive systems and are the true signature of a dispersive shock. The velocities of the two extrema of the oscillating region identify the so-called shock fan, which matches on both sides of the hydrodynamic (or dispersionless) limit of the original equation [5]. Dispersive shocks are ubiquitous phenomena that are found in a large variety of physical contexts, including Bose-Einstein condensates (BEC) [6][7][8][9][10][11][12][13][14][15], optics [16][17][18][19][20][21][22], dispersive gas dynamics [23], plasma physics [24,25] and granular systems [26].
One of the most complete frameworks for understanding DSWs is represented by the Whitham theory of modulation [27]. The main idea behind this approach is to identify different scales of evolution and to average-employing a suitable Whitham averaging procedure-over the fast oscillatory scale of the dynamics. Such an averaging technique yields a set of modulation equations that describe the spatio-temporal evolution of all the parameters ruling the DSW. In the Whitham approach, one of the most stringent conditions lies in the existence of a sufficiently large set of conservation laws, which should match the number of parameters necessary to represent the modulated periodic wave and hence the DSW. To date, the Whitham theory has been successfully applied to various integrable systems, which are characterized by an infinite set of invariants of motion and hence do not pose any limit on the applicability of this method. In particular, a complete application of the Whitham theory has been developed for the Korteweg de-Vries (KdV) [1,28] and for the defocusing nonlinear Schrödinger (NLS) [29,30] equations 3 (see also the papers by Lax and Levermore [31] and Jin et al [32], which address rigorously the small dispersion limit of the KdV and NLS equations, respectively). Importantly, however, the lack of sufficient conservation laws prevents the same approach from being applicable to a wide class of non-integrable systems. In this case, the application of Whitham averaging is still a topic of ongoing research. When the method is aimed at the limited goal of describing the key dynamical properties of DSWs, suitable changes of variables can be found which can lead to Whitham averaged equations (i.e. modulation equation) that can be effectively solved [33]. This approach uses a particular reformulation of the Gurevich-Pitaevski construction [1] and some general properties of averaged quantities in order to diagonalize the Whitham equations of motion at the edges of the shock fan region, thus providing analytic expressions for the edge velocities of the shock fan, the wavenumber and the peak amplitude of the DSW. It is worthwhile observing that, despite the same name, the non-integrable version of the Whitham theory is based on a completely different approach with respect to the classical formulation of this method for integrable systems. In the latter case, in fact, diagonal modulation equations follow from the Lax pair structure underlying the integrability, and are characterized by Riemann invariants belonging to the periodic solution of the initial nonlinear equation [34]. In the former approach, conversely, a completely different set of invariants is introduced, which rely on the soliton solutions of the initial model, thus yielding a completely diverse set of diagonal equations (see section 3). Such a non-integrable version of the Whitham theory has been applied to the study of DSW in photorefractive crystals [35].
Among several non-integrable systems of interest in physics, the cubic-quintic nonlinear Schrödinger (CQNLS) equation represents a universal model widely employed for studying nonlinear dynamics. Besides being regarded as a general model for superfluidity [36], it finds direct application in BEC and optics. In particular, in the latter context, the CQNLS model describes the paraxial evolution of a beam propagating in a centrosymmetric medium where the two lowest order nonlinear terms yield a refractive index change with the intensity |ψ| 2 , which follows the law n(|ψ| 2 ) = n 2 |ψ| 2 + n 4 |ψ| 4 [37,38], with the parameters n 2 and n 4 being experimentally measurable by two-wave coupling or Z-scan techniques [39,40]. While DSW occurs in defocusing media (n 2 < 0), the regime of major interest entails n 4 > 0, for which n(|ψ| 2 ) represents, e.g., the power expansion of any saturable cubic nonlinearity or the macroscopic contribution arising from cascading effects. On the other hand, in the field of BEC, the CQNLS (or Gross-Pitaevskii) equation accounts, in the framework of mean field theory, for elastic two-and three-body interactions [41]. While atom-atom (two-body) interactions are well understood quantitatively and demonstrated to be controllable to a large extent, the effect of higher-order atom interactions is the subject of current investigations [42,43], and reliable data on the quantitative impact of elastic three-body contributions are lacking.
While the study of DSW has been carried out with significant efforts in the case of the classical cubic NLS [18], the effect of quintic nonlinearities is fairly unexplored, while it can be important for designing new experiments. Recently, based on the hydrodynamic reduction of the CQNLS equation, we have shown that a rich DSW crossover scenario can occur as a result of the growing impact of quintic nonlinearities [44]. In this paper, we specifically focus on the Whitham approach in the case when the DSW is fully stable against all possible variations of the control parameters. We report a detailed derivation of the Whitham equations, and demonstrate through their solutions that the DSW is tunable to a large degree (in terms of shock fan and shock velocities) by acting on a simple parameter, namely the input optical intensity or the density in BEC. Quite remarkably, the tunability is accessible in the limit of a perturbative nonlinearity.

4
Theoretical predictions are then verified against numerical simulations, finding very good agreement with analytical estimates for the DSW evolution. The paper is organized as follows. Section 2 describes the model and recalls the properties of solutions of the dispersionless limit of the CQNLS equation. Section 3, conversely, encompasses the application of the Whitham theory to the analysis of DSW. Parameters for shock amplitude, velocity and wavenumber are given in sections 3.1 and 3.2. Finally, numerical results are presented and compared with theoretical predictions in section 4, which also discusses a possible experimental observation of this dynamics.

The model and dispersionless limit of the cubic-quintic nonlinear Schrödinger equation
We begin with the following adimensional CQNLS: with 1 and α being an external positive free parameter that defines the relative strength of the quintic nonlinear response. In optics, ψ represents a slowly varying optical beam envelope (|ψ| 2 is in units of the reference input intensity I 0 ), α = n 4 I 0 /n 2 , and = √ L nl /L d is fixed by the ratio between linear L d = k 0 nw 2 0 and nonlinear L nl = (k 0 n 2 I 0 ) −1 length scales, with w 0 being the characteristic width of the input wavepacket ψ, k 0 = 2π/λ is the wavenumber and n is the linear index. In BEC, ψ corresponds to the mean-field boson wavefunction in units of the square root of a reference density ρ ∞ , = 1/(ρ ∞ a 2 0 a) and α = ρ ∞ g 3 /g 2 , with g 2 and g 3 being two-and three-body nonlinear coefficients, respectively [12,45], a is the scattering length associated with two-atom interaction, a 0 = √h /(mω) is the harmonic oscillator characteristic width, m is the atomic mass and ω is the transverse trap frequency. By employing the Madelung transformation after straightforward algebra, we obtain the equivalent hydrodynamic formulation of equation (1): where The variables ρ and u correspond to the density and velocity, respectively, of an incompressible fluid following the Eulerian equations of motion (3) and (4), with equation (5) representing the associated equation of state. The dispersionless limit is obtained by imposing → 0 in equation (4): 5 We consider the Riemann problem associated with equation (6), i.e. the evolution of a step-like initial condition in the density, with This condition has the double advantage of being analytically tractable and experimentally accessible in both nonlinear optics or BEC, assuming that the intensity or density slowly decays to zero far from the region around the discontinuity that embraces the relevant dynamics, so as to keep finite the power or atom number. On the other hand, an abrupt variation of intensity or atom density can be easily impressed, in the latter case by exploiting, e.g. a strong focused beam impinging on the condensate [12,21]. In the general case of smooth inputs, DSWs form after a finite distance (or time in BEC) where the wave breaks exhibiting an infinite steepness (gradient catastrophe). In this respect the initial discontinuity in equation (8) can be thought of as a case of breaking at zero distance, which is sufficiently accurate for assessing how the quintic nonlinearity affects the parameters of the DSW. The case of a smooth input will be assessed numerically at the end of section 4.
Equation (6) can be solved by the introduction of the Riemann invariants: with which diagonalize the system (6), obtaining The quantities V ± = u ± ρ − 2αρ 2 represent the Riemann characteristic eigenvelocities. Whenever the constraint is fulfilled, the Riemann velocities are real, and equation (6) forms a hyperbolic system that admits analytical solutions [46]. Since equations (11) are invariant under the scale transformation x → C x and z → C z, with C being arbitrary constant, they can be expressed in terms of the self-similar variable ξ = x/z, obtaining System (13) admits the following two simple wave solutions, where one of the Riemann invariants is kept constant while the other one is free to vary: Equations (14) and (15) define implicitly the solution through the function ξ = V ± = f (ρ). A final step consists in matching equations (14) and (15) with the initial condition (8). Imposing ρ = ρ − 0 in the second of equation (15) and using the definition of λ + in equations (9), we find the relation which allows us to reexpress u as a function of ρ. Upon substitution of equation (16) into the first of equation (15), we obtain which is the complete solution for the left propagating wave, called the left simple wave (figure 1(a)). Imposing ρ = ρ − 0 in (17), we obtain the velocity of the point x e− (figure 1(a)), which defines the extremum of the left simple wave: For x > 0, conversely, we have ρ = ρ + 0 with u = 0. The solution that satisfies this condition is a right simple wave obtained by matching equation (14) with the input (equation (8)). This yields the right simple wave solution The velocity of the extremum point x e+ (see figure 1) of this right simple wave is found by imposing ρ = ρ + 0 in equation (19), obtaining The connection between the two simple waves is realized by introducing a constant density ρ i and velocity u i distribution (see figure 1). Values for ρ i and u i are found by imposing ρ = ρ i in equations (17)- (19) and u = u i in equation (19), thus obtaining The intermediate constant level ρ i , u i is bounded in x by two internal points x i± , which evolve with velocities ξ i± ( figure 1(a)). Figure 2 shows the behavior of the intermediate density versus α, as calculated by a numerical solution of equation (21). We observe that the density ρ i is almost constant when the control parameter α is varied, and exhibits only a modest decrease. Substituting equation (21) into the simple wave solutions (14) and (15), we obtain Figure 1 displays the behavior of the density ρ as follows from equations (17)- (19). As shown, the initial discontinuity decays into two distinct waves connected by an expanding region of constant density ρ i : the first is a rarefaction wave which moves toward the left (0 < x) and the second is a multivalued (broken) wave which moves toward the right (x > 0), owing to the fact that the point x i+ moves with velocity ξ i+ larger than the velocity ξ e+ of the point x e+ . The latter fact is indeed at the origin of the fact that a DSW can be formed in the region x e+ x x i+ where spontaneously forming oscillations regularize the overturn found here in the dispersionless limit. As seen from equations (18)-(22), a condition required for the formation of a multivalued region is ρ 0+ > 0. If this condition is not met, the velocity ξ e+ = 0 while ξ i+ < 0, which does not lead to any wave breaking for x 0 but conversely to the propagation of two simple waves for both x < 0 and x > 0.

Whitham analysis of the dispersive shock waves
We study the regularization of the multivalued region by employing the non-integrable form of the Whitham modulation theory, which is extensively discussed in [33]. Following this approach, the DSW is modeled as a periodic wave with slowly varying the amplitude and phase.
The velocity v of a generic point x(z) of the DSW can be represented as follows: 8 and depends on the following four variables: the amplitude a of the modulated wavetrain, the Whitham averaged density ρ and velocity u, and the wavenumber k = 2π/L associated with the spatial period L of the periodic wave. The latter defines in turn the Whitham average of any arbitrary quantity y [27]: Here Q(ρ) is a polynomial which arises when searching for periodic solutions of the starting model (1). Following the approach of [47], one can easily show that periodic waves traveling at velocity c, i.e. dependent on the retarded variable θ = x − cz, can be found for our starting model, as the Jacobian solutions of an ordinary differential equation (ODE) where Q(ρ) can be expressed in the general form with H being the Hamiltonian (total conserved energy) associated with the equivalent motion of an ideal particle of the coordinate ρ, and V (ρ) taking the form of a quartic (doublewell) potential. In general, two sets of solutions exist for a given value of energy H, each corresponding to the motion bounded within one of the two wells. When H approaches the value that pertains, for a given set of parameters, to the homoclinic separatrix of the system, the period L tends to infinity, and the two sets of periodic waves tend to a dark (a hole on a constant background) and an antidark (a bright disturbance on the same background) soliton, respectively, as extensively discussed in [47]. The Whitham theory that we discuss below does not necessitate writing explicitly such periodic solutions, whereas it makes use of the Whitham modulation equations in the specific limit for which the periodic wave of equation (1) reduces to the dark soliton. In this case it is convenient to parametrize the solution in terms of an additional parameter (besides the velocity c), namely the intensity (or density) ρ i , which corresponds to the soliton background in the limit of an infinite period. In terms of the two parameters ρ i , c, the potential V explicitly reads as (see [47] for details) where k 0 = c 2 2 + f (ρ i ). The region where the modulated periodic wave develops thus regularizing the multivalued structure discussed in the previous section-usually called the dispersive shock fan-is bounded by two different coordinates x e and x i , called the trailing edge and the leading edge of the fan, respectively. These edge coordinates move with velocities described by a reduced set of variables: Along the leading edge, in fact, the periodic wavetrain approaches a solitary wave, thus requiring k → 0, whereas along the trailing edge the periodic wave becomes a linear periodic solution with a → 0. The key point of the non-integrable formulation of the Whitham theory is the reduction of the Whitham equations to a suitable dispersionless-like form, which can be 9 effectively solved at the trailing/leading edge thanks to a reduced number of unknowns. As we will discuss in the next paragraphs, the solution of the Whitham equations at the shock fan boundaries will be sufficient to study the main properties of the DSW.

Dynamics at the trailing edge
The Whitham equations are obtained by averaging the conservation laws of the CQNLS equation. The general form of these equations can be represented as follows: Equations (29) and (30) represent the conservation of mass and momentum, respectively. A third equation is derived from the wavenumber conservation law [34, section 3.5]: with ω = vk. The system composed of equations (29)-(31) provides a set of three equations for four unknowns, and cannot be solved inside the shock fan region. However, equations (29)- (31) can be studied at the trailing/leading edge, where three unknowns are sufficient for describing the DSW dynamics. To this extent, we exploit the property [33] y(ρ, u)| a→0 = y(ρ, u), and take the limit a → 0 in equations (29) and (30). With the aid of equation (32), we obtain the following dispersionless-like form: Correspondingly, the wavenumber conservation law becomes where stands for the dispersion relation obtained from the linearization of equations (3)-(4) with ρ = ρ + ρ 1 exp i(kx − ωz) and u = u + u 1 exp i(kx − ωz), with |ρ 1 |, |u 1 | much smaller than their constant backgrounds ρ, u. Equations (33)-(36) form a closed system that can be integrated to derive an analytic description of the DSW at the trailing edge. Equations (33) and (34) are decoupled from equation (35), and can be independently solved for ρ(x, z), u(x, z). We also observe that equations (33) and (34) are identical to equations (3) and (4), except for the presence of averaged quantities instead of local ones. Therefore, all the results obtained in section 2 apply here. In particular, we have which is the counterpart of equation (19). Substituting equation (37) into equations (33)-(35), we obtain where and The system (38) admits two families of characteristic velocities: The family in the first of equation (42) does not depend on k and represents the equivalent of the first of equation (19), while the second of equation (42) determines the k-dependent velocity v 0 . By looking for solutions of the type k = k(ρ), equation (38) reduces to We then introduce a change of variable by defining By differentiating both terms of equation (44), with the aid of equation (41), we arrive at the ODE dγ (ρ) This equation must be solved in conjunction with the constraint k(ρ i ) = 0, which represents the condition for the degeneracy into a solitary wave at the leading edge. This yields the initial value for the ODE (45), which is expressed as γ (ρ i ) = 1. Once the solution γ (ρ) is found, we can calculate the wavenumber k + of the linear wavepacket as follows: as well as the velocity v 0 = dx e /dz of the DSW edge x e : We observe that equation (47) is not analytically solvable. However, a solution in the limit of a small α can be found by perturbative analysis. In particular, in the case α = 0 we have We then look for a perturbative solution of equation (45) expressed as γ =γ +γ 1 . After long but straightforward algebra, we find the first-order perturbative correction termγ 1 : We then expand equation (47) at first order in α and, by substituting the previous expression, we find that which provides an analytic solution to the problem.

Leading edge
At the leading edge, we need to solve equations (33)- (36) in the limit where the DSW approaches a soliton of the CQNLS, i.e. for k → 0. In this case, however, a direct generalization of the previous approach is not feasible, due to the singular limit of equation (35) when k tends to 0. In order to overcome the problem, we use the conjugate wavenumber transformation introduced in [33]. This approach employs a coordinate change in the system, defining a new amplitude:ã and a new wavenumberk = k/ã. The parameters (ρ, u,k,ã) provide a new set of variables for the Whitham system (33)- (36), where a new wave (called the conjugate wave) is defined throughk andã. In the leading edge, the DSW approaches a soliton solution of the CQNLS and, consequently, equation (26) develops a double root for ρ = ρ i [47]: where the roots ρ a , ρ m read as The next step is to obtain the relationship between the amplitude of the conjugate wave and the wavevector k of the original wave. To this extent, we use the exponential decay of intensity ρ i − ρ ∝ exp(−χ|x − v s z|), for |θ | = |x − v s z| → ∞, and expand Q near the double root ρ i , thus obtainingã In the non-conjugate representation, the term exp χ(x − v s z) is mapped into exp i(kx − χv s z); hence we have χ = k. This, in turn, leads toã = k, which states that at the leading edge the new 12 amplitudeã coincides with the wavenumber k of the original (i.e. non-conjugated) wave. In the representation withã andk, therefore, the condition k → 0 is equivalent toã → 0, which does not generate any singularity in equation (35). We can now repeat the same analysis developed in section 3.1; in particular, the dispersion relation for the conjugated system reads as Substituting k =kã into (40), we obtaiñ which is satisfied for k = 0 (i.e.ã = 0). When k → 0, we can expand the termk(ã z +ω x ) for a → 0 and obtain the wavenumber conservation law at the leading edge: with˜ (ρ, k) that reads as Equation (59) is obtained by substituting equation (37) into equation (56). We now proceed to the solution of equation (33), following the same type of analysis discussed in section 3.1. In particular, by employing equation (37), we obtain with V (ρ) defined from equation (40). From equations (58)-(60), we obtain the two families of characteristic velocities: The first of equations (61) represents the simple wave solution of the Whitham system, while the second is the equivalent of equation (28). By seeking the solution of the type k = k(ρ) and by employing equation (56), we obtain Then by defining the quantity we obtain the analogue of equation (45) in the variableγ (ρ): This equation must be solved with the initial condition k(ρ + 0 ) = 0, which implies thatγ (ρ + 0 ) = 1. Once the solution of equation (64) is known, we can determine the velocity of the edge x i+ from equations (61) and (59): 13 The amplitude a, conversely, is calculated from the density jump Analytic solutions are found by perturbative analysis, as developed in section 3.1. In particular, we begin by looking for solutions of the type γ =γ +γ 1 withγ given by equation (48) and We then perturbatively expand equation (65) and find that

Numerical experiments
In order to check the outcome of our analysis, we have performed numerical simulations by developing a fourth-order, finite difference integration scheme for equation (1). At the input ψ(x, z = 0) = ψ 0 , we considered a sharp step-like variation generated by one edge of a super-Gaussian pulse distribution ψ 0 = 1 + A exp[−(x/w) p ], with amplitude A = ρ − 0 − 1, waist w = 1.5 and exponent p = 60. We have then monitored the pulse propagation for different values of the parameter α. Figure 3 summarizes the results of our numerical analysis, showing also a comparison with the analytical prediction arising from the Whitham theory. At α = 0.01 (figure 3(a)), the system develops a single DSW for x > 0 and a simple rarefaction wave propagating at x < 0, in agreement with the solution of the Riemann problem of the dispersionless limit (see section 2). Figure 3 also reports the velocities of the shock fan v s and v 0 (solid lines) calculated from equations (47) and (65), as well as the velocities ξ i− , ξ e− , ξ i+ and ξ e+ (dashed lines) of the simple waves computed from equations (18)- (20) and (22). While for the rarefaction waves (x < 0) the dispersionless limit provides excellent agreement with numerics (see figures 3(a) and (b)), it significantly fails to predict the system dynamics in the shock region x > 0 (see figures 3(a) and (b), dotted lines for x > 0). In the latter case, in fact, regularization of the wave breaking is mainly achieved through dispersion-which is neglected in the hydrodynamic approximation-and can thus be described only within the Whitham approach that accounts for all the dynamic effects of the starting CQNLS equation (1). As shown in figures 3(a) and (b), the predictions arising from the Whitham theory are extremely accurate when compared to numerical results. When α is increased in the model, we observe a steering of the main direction of propagation of the shock fan toward the center x = 0 (compare figure 3(a) with figure 3(b)), characterized by unequal changes of the velocities v s and v 0 . This can be analytically predicted on the basis of equations (50) and (68), which both show a decreasing of the shock wave velocity for increasing α. In particular, the velocity v s changes more smoothly with α when compared with v 0 . Therefore, this generates a compression of the dispersive shock toward x = 0. A snapshot of the DSW at z = 10 for α = 0.1 is reported in figure 3(c), which shows the formation of the fast modulated wavetrain that regularizes the discontinuity for x > 0. The modulated region, in particular, starts from a linear background at x e+ and matches on the other edge v s a dark soliton solution of the CQNLS (which is always stable for any value of the parameters [47]). Further increase of the propagation distance z does not qualitatively alter this scenario, but simply increases the number of oscillations inside the shock fan due to  (66)) as a function of α, overlapped with the numerical results. Here the initial step-like is such that ρ − 0 = 1.5 and ρ + 0 = 1.
an augmented spatial extension of the DSW along the x-axis. Figure 3(c) clearly shows the nonlinear steering of the shock toward the origin, as soon as the strength of the coefficient of the quintic terms α increases. The latter, in turn, depends on the field intensity I 0 (see section 2) and allows complete controllability of the DSW by acting on the field power P = I 0 dxρ. It is worthwhile remarking the small nature of the variation α ≈ 0.1 considered, equivalent to an intensity increase I 0 of ≈ 10%, which is already sufficient to yield a significant steering in the direction of propagation of the DSW. Figure 3(d) finally compares the behavior of the dark soliton amplitude a, as calculated from equation (66), with numerical estimates obtained by simulations. Such an analysis shows that the amplitude of the soliton is particularly stable, and does not present a significant variation with α, while matching very well with the theoretical predictions.
We complete our numerical analysis by studying the dynamics of an initial Gaussian pulse over a pedestal, by taking ψ 0 = 1 + exp (−x 2 ). This constitutes the simplest accessible input condition for experiments in both nonlinear optics and BEC systems. Figure 4 summarizes our results for α = 0.01 and 0.1. The system develops two symmetrically displaced DSWs which develop at finite distance from symmetrically located points of breaking along the tails of the Gaussian pulse. In the proximity of x = 0, conversely, the dynamics shows the propagation of two simple waves, such as in the previous case, which originate a plateau in the field intensity ρ. In the same fashion as in the sharp discontinuous case, an increasing of α is accompanied by a steering of the DSW toward the origin x = 0, as well as by a shock compression due to the different variation of velocity over the fan region (see figures 4(c) and (d)). These simulations show that the DSW dynamics is not qualitatively altered by the presence of a smooth input function, and is therefore a very robust phenomenon supported by the presence of competing nonlinearities in the system.

Conclusion
In conclusion, we have presented a theoretical analysis of the dispersive shocks arising in a nonlinear system encompassing high-order interactions of opposite signs. By employing a rigorous approach based on a non-integrable formulation of the Whitham theory of modulation, we developed a fully analytic study of the DSW dynamics. Even in the perturbative regime, the system shows a rich physical behavior, characterized by the formation of a DSW with completely controllable properties, characterized by a gradual steering and compression of the shock fan. Theoretical predictions are confirmed by numerical results. Finally, simulations carried out with input profiles more easily achievable in typical nonlinear optics and BEC experiments demonstrate the same qualitative behavior and highlight the robustness of this dynamics. This work is expected to stimulate further theory and experiments in both nonlinear optics and superfluidity, where the CQNLS finds different important applications.