11institutetext: Artificial Intelligence & Analytics, Data Science &
Artificial Intelligence, R&D, AstraZeneca, Gaithersburg, Maryland, USA
11email<EMAIL_ADDRESS>
# Personalized Step Counting Using Wearable Sensors: A Domain Adapted LSTM
Network Approach

Arvind Pillai  Halsey Lea  Faisal Khan  Glynn Dennis

###### Abstract

Activity monitors are widely used to measure various physical activities (PA)
as an indicator of general health. Similarly, real-time monitoring of
longitudinal trends in step count has significant clinical potential as a
personalized measure of disease related changes in daily activity. However,
inconsistent step count accuracy across vendors, body locations and individual
gait differences limits clinical decision making. The tri-axial accelerometer
inside PA monitors can be exploited to improve step count accuracy across
devices and individuals. In this study, we hypothesize: (1) raw tri-axial
sensor data can be modelled to create reliable and accurate step count, and
(2) a generalized step count model can then be efficiently adapted to each
person’s unique gait pattern using very little new data. Firstly, open-source
raw sensor data was used to construct a long short term memory (LSTM) deep
neural network to predict step count accurately ($\sim$98%). Then we generated
a new, fully independent dataset using a different device and different
subjects. Finally, a small amount of subject-specific data was domain adapted
to produce personalized models with high individualized step count accuracy
($\sim$98%). These results suggest models trained using large freely available
datasets can be fine-tuned on patient populations lacking large historical
data sets.

###### Keywords:

Step count Accelerometer Deep Learning.

##  1 Introduction

In 2008 the US government released its first-ever evidence-based guidance on
PA, which acknowledged very strong evidence that physically active people have
higher levels of health-related fitness, a lower risk profile for developing
cardiovascular diseases, and lower rates of various chronic diseases than
people who are not active [1]. A simple and inexpensive option for monitoring
steps and cardiovascular health is the physical activity monitor. Step count
is a quantifiable metric which demonstrates individuals taking $\geq$5,000
steps/day had substantially lower prevalence of adverse cardio-metabolic
health indicators than those taking lower amounts [2]. In addition, step count
has been shown to be inversely related to risk of progression to diabetes [3].
Collectively, reduced physical activity as measured by step count is a
consistent predictor of death in chronic heart failure, possibly surpassing
traditional laboratory-based exercise tests.

Recent advances in deep learning have enabled more sophisticated modelling of
patterns in sensor data to classify human physical activities. LSTM based
architectures have shown their efficacy in accurately predicting step count
[4, 5]. Herein we describe a personalized step count algorithm that adapts to
a subject’s unique gait pattern with very little new subject-specific data
collected from a different wearable device. We consider an LSTM deep learning
algorithm to model human steps from two different wrist-worn devices
containing accelerometers and gyroscopes. We observed a general LSTM model
trained on publicly available data performs well for most subjects, but
exhibited lower step count accuracy for some individuals. We hypothesized that
the observed decrease in step count accuracy may be due to individual subject
gait differences and change in devices. Therefore, we show that an adaptive
step count algorithm can be used to train personalized models to obtain
accurate step count across devices and subjects.

##  2 Materials and Methods

In this paper, we use two datasets: (1) the ”Pedometer” dataset [6], collected
at Clemson University using a Shimmer3 device from 30 subjects and (2) a fully
independent, analogous ”AZ” dataset, collected at AstraZeneca using an
ActiGraph GT9X Link from 5 different subjects. The participants were asked to
perform the 6-minute walk test in an indoor looped setting. Next, the steps
and distance were manually annotated using video recordings from an Iphone XR.
Morever, the rationale for developing this dataset was to test the ability of
an LSTM step count algorithm to learn features that generalize across
different device manufacturers and subjects. Accelerometer and gyroscope data
from the wrist worn devices were used to model step count as a classification
problem with left and right step as the two classes. The sensor data was
resampled to 15Hz and windowed into $0.4\overline{6}$s readings based on
tuning experiments and literature [7].

The deep learning model architecture consists of a many-to-one 256 unit LSTM
followed by two fully connected layers of 512 and 256 units, respectively,
with dropouts and activated using ReLU. Finally, a softmax layer was used to
predict left or right step, and every transition between the two classes is
counted as a step. The input features to the deep learning model architecture
consists of individual windows with 3 readings from a tri-axial
accelerometer(gs) and 3 from a tri-axial gyroscope(deg/sec). In order to
personalize our model without the need to train from scratch, we investigated
different variations of Transfer Learning. Domain adaptation is one such
paradigm where the source and target probability distributions are different,
i.e. given a source feature space $X_{S}$ and target feature space $X_{T}$,
source label space $Y_{S}$ and target label space $Y_{T}$, we have
$X_{S}=X_{T}$, $Y_{S}=Y_{T}$, and $P(X_{S})\neq P(X_{T})$ [8]. In the context
of our problem, $X_{G}=X_{P}^{S}$, $Y_{G}=Y_{P}^{S}$, and $P(X_{G})\neq
P(X_{P}^{S})$, where G represents a general model trained using Pedometer
dataset (source) and P represents a personalized model for every subject S in
the AZ dataset (target).

The metrics chosen to evaluate our model are class accuracy and step count as
shown in (1) and (2).

$accuracy_{class}=\left(\frac{N_{R}+N_{L}}{N_{T}}\right)*100$ (1)
$accuracy_{steps}=\left(1-\left|\frac{step\\_count_{predicted}-step\\_count_{groundtruth}}{step\\_count_{groundtruth}}\right|\right)*100$
(2)

where $N_{R}$ and $N_{L}$ are number of correctly classified right and left
step windows respectively, and $N_{T}$ is the total number of steps.
$step\\_count_{groundtruth}$ is obtain using visually confirmed annotation
from the datasets and $step\\_count_{predicted}$ is the step count predicted
by the model.

##  3 Results

To establish a baseline, we obtained data from Pedometer dataset [6] and
implemented Piece-wise Aggregate Approximation (PAA) + merge approach
described in [9].

After establishing a baseline, the LSTM deep learning model was designed as
described in the previous section. There are three major experimental setups
to assess the performance of our model. Firstly, a general model was trained
and tested using only the Pedometer dataset (Shimmer3) using Leave-two-
subject-out cross-validation (CV), we call this LSTM-General I. Next, this
pre-trained general model was validated using the entire AZ dataset (ActiGraph
GT9X Link) as a test set, which is documented as experiment LSTM-General II.
Finally, a personalized model, LSTM-Personalized was domain adapted using a
2-step training process from the pre-trained general model using only 30
seconds of subject-specific data, and the test dataset was created by removing
this information from each subject. As shown in Table 1, the median
$accuracy_{class}$ and median $accuracy_{steps}$ among subjects were recorded
for the three settings. Since, LSTM-General I generates different models for
each CV fold, we calculate the mean of medians and use the best pre-trained
model.

The model was tuned based on window size, sensor modalities, device location,
and number of units in each layer. In addition to training with different
combination of layers, we also examined different LSTM architectures including
stacked, bidirectional, and stateful.

Table 1: Summary of $\mathbf{accuracy_{class}}$ and $\mathbf{accuracy_{steps}}$ for all experimental setups. Approach | Training Data | Testing Data | $\mathbf{accuracy_{class}}$ | $\mathbf{accuracy_{steps}}$
---|---|---|---|---
PAA + Merge (Baseline) | N/A | Pedometer | N/A | 89.73
PAA + Merge (Baseline) | N/A | AZ | N/A | 86.05
LSTM-General I | Pedometer | Pedometer | 85.98 $\pm\,2.46$ | 98.63 $\pm\,2.37$
LSTM-General II | Pedometer | AZ | 60.01 | 94.49
LSTM-Personalized | Pedometer + $\sim$AZ | AZ | 91.76 | 98.81

As shown in Table 1, LSTM-General I improves $accuracy_{steps}$ by $\sim$10%
from the baseline, and the standard deviation indicates a left-skewed
distribution. Moreover, the baseline tested is an approach that does not
classify steps, so $accuracy_{class}$ metric is not applicable. Additionally,
there is a decrease in performance of LSTM-General II which is indicated by a
25.97% and 4.14% drop in $accuracy_{class}$ and $accuracy_{steps}$
respectively, which shows lack of robustness in handling different devices.
Next, we observe LSTM-Personalized outperforms LSTM-General II as indicated by
a 31.75% and 4.32% increase in $accuracy_{class}$ and $accuracy_{steps}$
respectively. This implies only a small amount of labelled data (represented
as $\sim$AZ in Table 1) is necessary to improve performance.

(a) (b)

Figure 1: Comparison of experiments at subject level: (a) $accuracy_{class}$
box plot. (b) $accuracy_{steps}$ box plot.

.

From the boxplots shown in Fig. 1., drastic decreases in performance of LSTM-
General II in case of outliers is also adequately handled by LSTM-
Personalized. Additionally, the median performance of LSTM-Personalized and
LSTM-General I are similar in terms of both the metrics, which suggests
information from different devices and different subjects can be handled by
using personalized models.

##  4 Conclusion

In this paper, we proposed a personalized step counting algorithm built using
an LSTM network and tuned with subject-specific data through domain
adaptation. This shows that data collected from different subjects and devices
can be modelled using our algorithm. Additionally, results suggest very little
labelled data is necessary to construct personalized models. This work
suggests deep transfer learning has significant potential to achieve
clinically acceptable accuracy in the assessment of physical activity for
multiple disease scenarios. Moreover, accurately identifying class labels is
important as it helps in tasks like distance estimation, indoor localization,
and human activity recognition. Therefore, in the future, we plan to conduct a
large-scale study in patient populations to evaluate other wearable sensor
related tasks, disease endpoints, and unsupervised domain adaptation.

## References

  * [1] Tudor-Locke C: Steps to Better Cardiovascular Health: How Many Steps Does It Take to Achieve Good Health and How Confident Are We in This Number? Curr Cardiovasc Risk Rep. 4(4) (2010)
  * [2] Schmidt MD, Cleland VJ, Shaw K, Dwyer T, Venn AJ: Cardiometabolic risk in younger and older adults across an index of ambulatory activity. Am J Prev Med. 37(4) 278–284 (2019)
  * [3] Kraus WE, Yates T, Tuomilehto J, et al: Relationship between baseline physical activity assessed by pedometer count and new-onset diabetes in the NAVIGATOR trial. BMJ Open Diabetes Research and Care. 6e000523 (2018)
  * [4] M. Edel and E. Köppe: An advanced method for pedestrian dead reckoning using BLSTM-RNNs. International Conference on Indoor Positioning and Indoor Navigation (IPIN). 1–6 (2015)
  * [5] Ziyi C: An LSTM Recurrent Network for Step Counting. arXiv preprint arXiv:1802.03486 (2018)
  * [6] R. Mattfeld, E. Jesch and A. Hoover: A new dataset for evaluating pedometer performance. 2017 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), Kansas City, MO. 865–869 (2017)
  * [7] Pachi, A. & Ji, Tianjian.s: Frequency and velocity of people walking. Structural Engineer 83, 36–40 (2015)
  * [8] Bleiweiss, A.: LSTM Neural Networks for Transfer Learning in Online Moderation of Abuse Context. ICAART 2, 112-122 (2019)
  * [9] Yuanyuan F et al.: Comparison of tri-axial accelerometers step-count accuracy in slow walking condition. Gait & Posture. 11–16 (2017)


* Villamizar et al. (2022) Hugo Villamizar, Marcos Kalinowski, and Helio Lopes. 2022\.  A Catalogue of Concerns for Specifying Machine Learning-Enabled Systems.  (2022).  http://arxiv.org/abs/2204.07662
  * Wachter et al. ([n. d.]) Sandra Wachter, Brent Mittelstadt, and Chris Russell. [n. d.].  Bias Preservation in Machine Learning : The Legality of Fairness Metrics Under EU Non- Discrimination Law.  _West Virginia Law Review, Forthcoming_ ([n. d.]), 1–51.  https://ssrn.com/abstract=3792772
  * Wang et al. (2017) Yuxuan Wang, Pascal Getreuer, Thad Hughes, Richard F. Lyon, and Rif A. Saurous. 2017.  Trainable frontend for robust and far-field keyword spotting.  _ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings_ 1 (2017), 5670–5674.  https://doi.org/10.1109/ICASSP.2017.7953242
  * Warden (2018) Pete Warden. 2018\.  Speech Commands: A Dataset for Limited-Vocabulary Speech Recognition. (2018).  https://arxiv.org/abs/1804.03209
  * Yang et al. (2018) Tien Ju Yang, Andrew Howard, Bo Chen, Xiao Zhang, Alec Go, Mark Sandler, Vivienne Sze, and Hartwig Adam. 2018\.  NetAdapt: Platform-aware neural network adaptation for mobile applications.  _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ 11214 LNCS (2018), 289–304.  https://doi.org/10.1007/978-3-030-01249-6{_}18
  * Zhang et al. (2017) Yundong Zhang, Naveen Suda, Liangzhen Lai, and Vikas Chandra. 2017.  Hello edge: Keyword spotting on microcontrollers.  _arXiv_ (2017), 1–14.

##  Appendix A Appendix

###  A.1. Experiment Setup: Datasets

Google Speech Commands keyword classes: ’bed’:0, ’bird’:1, ’cat’:2, ’dog’:3,
’down’:4, ’eight’:5, ’five’:6, ’four’:7, ’go’:8, ’happy’:9, ’house’:10,
’left’:11, ’marvin’:12, ’nine’:13, ’no’:14, ’off’:15, ’on’:16, ’one’:17,
’right’:18, ’seven’:19, ’sheila’:20, ’six’:21, ’learn’:22, ’stop’:23,
’three’:24, ’tree’:25, ’two’:26, ’up’:27, ’wow’:28, ’yes’:29, ’zero’:30,
’backward’:31, ’follow’:32, ’forward’:33, ’visual’:34

MSWC English keyword classes: ’about’: 0, ’after’: 1, ’also’: 2, ’been’: 3,
’could’: 4, ’first’: 5, ’from’: 6, ’have’: 7, ’however’: 8, ’just’: 9, ’know’:
10, ’like’: 11, ’many’: 12, ’more’: 13, ’most’: 14, ’only’: 15, ’other’: 16,
’over’: 17, ’people’: 18, ’said’: 19, ’school’: 20, ’some’: 21, ’that’: 22,
’they’: 23, ’this’: 24, ’three’: 25, ’time’: 26, ’used’: 27, ’were’: 28,
’what’: 29, ’when’: 30, ’will’: 31, ’with’: 32, ’would’: 33, ’your’: 34

MSWC German keyword classes: ’aber’: 0, ’alle’: 1, ’auch’: 2, ’dann’: 3,
’dass’: 4, ’diese’: 5, ’doch’: 6, ’durch’: 7, ’eine’: 8, ’gibt’: 9, ’haben’:
10, ’hauptstadt’: 11, ’heute’: 12, ’hier’: 13, ’immer’: 14, ’jetzt’: 15,
’kann’: 16, ’können’: 17, ’mehr’: 18, ’muss’: 19, ’nach’: 20, ’nicht’: 21,
’noch’: 22, ’oder’: 23, ’schon’: 24, ’sein’: 25, ’sich’: 26, ’sind’: 27,
’wenn’: 28, ’werden’: 29, ’wieder’: 30, ’wird’: 31, ’wurde’: 32, ’zwei’: 33,
’über’: 34

MSWC French keyword classes: ’alors’: 0, ’aussi’: 1, ’avec’: 2, ’bien’: 3,
’cent’: 4, ’cette’: 5, ’comme’: 6, ’c’est’: 7, ’dans’: 8, ’deux’: 9, ’donc’:
10, ’elle’: 11, ’fait’: 12, ’huit’: 13, ’mais’: 14, ’mille’: 15, ’monsieur’:
16, ’même’: 17, ’nous’: 18, ’numéro’: 19, ’plus’: 20, ’pour’: 21, ’quatre’:
22, ’saint’: 23, ’sept’: 24, ’soixante’: 25, ’sont’: 26, ’tout’: 27, ’trois’:
28, ’très’: 29, ’vingt’: 30, ’vous’: 31, ’également’: 32, ’était’: 33, ’être’:
34

MSWC Kinyarwanda keyword classes: ’abantu’: 0, ’ariko’: 1, ’avuga’: 2, ’bari’:
3, ’benshi’: 4, ’buryo’: 5, ’cyane’: 6, ’gihe’: 7, ’gukora’: 8, ’gusa’: 9,
’hari’: 10, ’ibyo’: 11, ’icyo’: 12, ’igihe’: 13, ’imana’: 14, ’imbere’: 15,
’kandi’: 16, ’kuba’: 17, ’kugira’: 18, ’kuko’: 19, ’kuri’: 20, ’mbere’: 21,
’muri’: 22, ’ndetse’: 23, ’neza’: 24, ’ntabwo’: 25, ’nyuma’: 26, ’perezida’:
27, ’rwanda’: 28, ’ubwo’: 29, ’umuntu’: 30, ’umwe’: 31, ’yagize’: 32, ’yari’:
33, ’yavuze’: 34

###  A.2. Statistical Analysis: Design Choices Arising During Model Training

Model 1: First factorial ANOVA interaction model for model training design
choices

⬇

model_inital = ’metric ~ C(dataset_name, Sum)+C(model_arch,
Sum)C(resample_rate, Sum)+C(mfccs, Sum)+C(mel_bins, Sum)+C(frame_length,
Sum)+C(frame_step, Sum)+C(window_fn, Sum)+C(dataset_name, Sum)*C(model_arch,
Sum)*C(resample_rate, Sum)*C(mfccs, Sum)*C(mel_bins, Sum)+C(dataset_name,
Sum)*C(model_arch, Sum)*C(resample_rate, Sum)*C(frame_length,
Sum)*C(frame_step, Sum)*C(window_fn, Sum)’

Model 2: Final factorial ANOVA interaction model for effect of model training
design choices on MCC

⬇

model_final_mcc = ’mcc ~ C(dataset_name, Sum)+C(model_arch,
Sum)+C(resample_rate, Sum)+C(mfccs, Sum)+C(mel_bins, Sum)+C(dataset_name,
Sum)*C(resample_rate, Sum)*C(mfccs, Sum)+C(model_arch, Sum)*C(mfccs,
Sum)*C(mel_bins, Sum)+C(dataset_name, Sum)*C(model_arch, Sum)*C(mfccs,
Sum)+C(frame_length, Sum)+C(frame_step, Sum)+C(model_arch,
Sum)*C(frame_length, Sum)*C(frame_step, Sum)’

Model 3: Final factorial ANOVA interaction model for effect of model training
design choices on reliability bias

⬇

model_final_bias = ’reliability_bias ~ C(dataset_name, Sum)+C(model_arch,
Sum)+C(resample_rate, Sum)+C(dataset_name, Sum)*C(resample_rate, Sum)+C(mfccs,
Sum)+C(mel_bins, Sum)+C(dataset_name, Sum)*C(model_arch, Sum)*C(mfccs,
Sum)+C(dataset_name, Sum)*C(mel_bins, Sum)+C(frame_length, Sum)’

###  A.3. Impact of Pre-processing parameters

Figure 13. Effect of MFCC dimensions on accuracy and reliability bias for
16kHz CNN models. Models without MFCC features, i.e. models that directly use
log Mel spectrograms as input features, tend to perform worse than those that
use MFCC features. Figure 14. Effect of log Mel spectrogram dimensions (# Mel
fbanks) on accuracy and reliability bias, disaggregated by input feature type
for CNN architectures. The number of Mel filter banks clearly impacts models
that directly use log Mel spectrograms as input features.

###  A.4. Statistical Analysis: Design Choices Arising During Model
Optimization

Model 4: First factorial ANOVA interaction model for pruning hyper-parameters

⬇

model_inital = ’delta_metric ~ mcc_baseline + reliability_bias_baseline +
C(dataset_name, Sum)+C(model_arch, Sum)+C(resample_rate,
Sum)+C(pruning_learning_rate, Sum)+C(pruning_schedule,
Sum)+C(pruning_frequency, Sum)+C(pruning_final_sparsity, Sum)+C(dataset_name,
Sum)*C(model_arch, Sum)*C(resample_rate, Sum)*C(pruning_learning_rate,
Sum)*C(pruning_schedule, Sum)*C(pruning_frequency,
Sum)*C(pruning_final_sparsity, Sum)’

Model 5: Final factorial ANOVA interaction model for effect of pruning design
choices on change in MCC

⬇

model_final_delta_mcc = ’delta_mcc ~ mcc_baseline + reliability_bias_baseline
+ C(dataset_name, Sum)+C(model_arch, Sum)+C(resample_rate,
Sum)+C(pruning_learning_rate, Sum)+C(pruning_schedule,
Sum)+C(pruning_frequency, Sum)+C(pruning_final_sparsity, Sum)+C(model_arch,
Sum)*C(pruning_learning_rate, Sum)*C(pruning_schedule,
Sum)*C(pruning_final_sparsity, Sum)+C(dataset_name, Sum)*C(model_arch,
Sum)*C(resample_rate, Sum)*C(pruning_final_sparsity, Sum)+C(dataset_name,
Sum)*C(pruning_learning_rate, Sum)*C(pruning_schedule, Sum)+C(dataset_name,
Sum)*C(pruning_schedule, Sum)*C(pruning_final_sparsity, Sum)+C(dataset_name,
Sum)*C(pruning_learning_rate, Sum)*C(pruning_final_sparsity,
Sum)+C(resample_rate, Sum)*C(pruning_learning_rate,
Sum)*C(pruning_final_sparsity, Sum)’

Model 6: Final factorial ANOVA interaction model for effect of pruning design
choices on change in reliability bias

⬇

model_final_delta_bias = ’delta_reliability_bias ~
mcc_baseline+reliability_bias_baseline+C(dataset_name, Sum)+C(model_arch,
Sum)+C(resample_rate, Sum)+C(pruning_learning_rate, Sum)+C(pruning_schedule,
Sum)+C(pruning_frequency, Sum)+C(pruning_final_sparsity, Sum)+C(dataset_name,
Sum)*C(model_arch, Sum)*C(resample_rate, Sum)*C(pruning_learning_rate,
Sum)+C(dataset_name, Sum)*C(pruning_learning_rate,
Sum)*C(pruning_final_sparsity, Sum)+C(pruning_schedule,
Sum)*C(pruning_final_sparsity, Sum)+C(model_arch,
Sum)*C(pruning_final_sparsity, Sum)’

###  A.5. Model Selection After Pruning

Figure 15. MCC (accuracy) score and reliability bias of models after pruning.
Models were selected based on three selection criteria: high accuracy, low
bias, and high accuracy + low bias. After pruning MCC is greatest for models
selected with a criteria that considers high accuracy. Similarly, reliability
bias is lower for criteria that consider low bias.


# Modulation theory for soliton resonance and Mach reflection

Samuel J. Ryskamp  Department of Applied Mathematics, University of Colorado,
Boulder, CO<EMAIL_ADDRESS>Mark A. Hoefer  Department of
Applied Mathematics, University of Colorado, Boulder, CO Gino Biondini
Department of Mathematics and Department of Physics, State University of New
York, Buffalo, NY

###### Abstract

Resonant Y-shaped soliton solutions to the Kadomtsev-Petviashvili II (KPII)
equation are modelled as shock solutions to an infinite family of modulation
conservation laws. The fully two-dimensional soliton modulation equations,
valid in the zero dispersion limit of the KPII equation, are demonstrated to
reduce to a one-dimensional system. In this same limit, the rapid transition
from the larger Y soliton stem to the two smaller legs limits to a travelling
discontinuity. This discontinuity is a multivalued, weak solution satisfying
modified Rankine-Hugoniot jump conditions for the one-dimensional modulation
equations. These results are applied to analytically describe the dynamics of
the Mach reflection problem, V-shaped initial conditions that correspond to a
soliton incident upon an inward oblique corner. Modulation theory results show
excellent agreement with direct KPII numerical simulation.

##  1 Introduction

Mach reflection occurs when a sufficiently large amplitude line soliton or
classical shock interacts with a barrier at a sufficiently small angle. A
Y-shaped triad is formed consisting of two smaller amplitude solitons or
shocks and a larger ”Mach” stem perpendicular to the barrier (see figure 1).
This phenomenon was first reported experimentally in the context of shallow
water solitons impinging on a corner in J. Scott Russell’s seminal 1845 paper
1845 [1], while later in 1875 Ernst Mach observed his eponymous phenomenon
arising through interacting shocks in gas dynamics [2, 3]. In this paper, we
present a method for modelling Y-shaped and X-shaped solitons as shock
solutions to soliton modulation equations.

Line solitons are special traveling wave solutions of the Kadomstev-
Petviashvili (KP) equation

$(u_{t}+uu_{x}+\epsilon^{2}u_{xxx})_{x}+u_{yy}=0,$ (1)

where $\epsilon>0$ is a parameter representing the strength of dispersion.
This version (1) with $+u_{yy}$ is also known as KPII [4], a completely
integrable equation [5] that arises in surface water waves [6, 7], internal
water waves [8, 9], and ion-acoustic waves in plasma [10]. It admits a two-
parameter family of stable line soliton solutions,

$u(x,y,t)=a\operatorname{sech}^{2}\left(\sqrt{\frac{a}{12}}\frac{x+qy-
ct}{\epsilon}\right),\qquad c=\frac{a}{3}+q^{2},$ (2)

where $a$ is the soliton peak amplitude, $q=\tan{\varphi}$, $\varphi$ being
the counterclockwise angle between the soliton and the $y$-axis (see figure
2$(a)$), and $\epsilon$ is proportional to the soliton’s width in the
direction of propagation. The KP equation (1) is a two-dimensional
generalization of the Korteweg-deVries (KdV) equation and therefore a natural
framework for the study of line solitons; when $q=0$ in (2), the well-known
KdV soliton is recovered.

The KP equation (1) also admits a three-parameter family of Y soliton
solutions that consist of a resonant triad of line solitons stably joined
together (see figure 2$(b)$) and represent the long-time asymptotic state for
line solitons after Mach reflection. Objects of much study, Y-shaped
configurations have been observed in numerics [11, 12, 13, 14, 15, 16, 17, 18,
19], physical experiments [1, 20, 21, 22, 23, 16, 24, 25, 26], and field
observations [6, 27, 28]. KP Y solitons were first described mathematically by
John Miles [29, 30], whose results were generalized by the discovery of more
complex KP exact solutions that contain Y-shaped substructures [31, 32, 33,
14, 34, 35]. Numerical studies have also identified approximately resonant
Y-shapes arising in other similar equations [36, 37, 38, 39].

Figure 1: $(a)$ Mach reflection of a soliton or shock impinging on a
compressive corner with sufficiently small angle $\varphi$. $(b)$ Mach
reflection also occurs for solitons with an initial inward V-shape, which
under the KP equation (1) is identical to the left figure for the purposes of
analysis. Figure 2: $(a)$ A line soliton exact solution (2) to the KP equation
(1), where $a$ represents the max amplitude, $q=\tan{\varphi}$, and the width
is proportional to $\epsilon$. $(b)$ A resonant Y soliton solution to the KP
equation (1) with physical variables $(a_{i},q_{i})$ of the $i^{\rm th}$
soliton labeled. The physical variables are determined by three phase
variables according to the conditions (6). This Y soliton corresponds to the
local behaviour of the bottom half of figure 1$(b)$.

In this paper, we study the Y soliton in the zero dispersion (i.e.
$\epsilon\to 0$) limit of the KP equation (1). In this limit, which is
equivalent to the limit of long time $t$ and large spatial $(x,y)$ scales for
some kinds of data, line solitons (2) vanish in width to become localized to a
line in the $x$-$y$ plane. When a line soliton is modulated, its parameters
$(a,q)$ in the zero dispersion limit evolve according to Whitham modulation
equations [40, 41, 42, 43, 44]. Recent advances in Whitham modulation theory
for the KP equation (1) [45, 46, 47] have proven to be effective for studying
the asymptotic evolution of line solitons interacting due to various initial
conditions [48, 49, 50], some of the first applications of Whitham theory to a
governing equation in two spatial and one temporal dimension.

In this work, we obtain two main results by examining the zero dispersion
limit of line solitons and Y solitons. The first result is that, in this
limit, all continuous line soliton modulations become essentially one-
dimensional in the $x$-$y$ plane. We demonstrate that $a$ and $q$ are
constrained so that if one knows the parameters’ evolution in one spatial
dimension, their evolution in the other dimension can be determined. This
result will allow us to reduce the (2+1)-dimensional soliton modulation
equations to a (1+1)-dimensional system.

The second and primary result of this paper is that the KP Y soliton can be
interpreted as a shock solution to the soliton modulation equations. In the
zero dispersion limit, the three-way intersection of line solitons in the
centre of the Y soliton reduces to a single point of discontinuity in the
modulation parameters. The velocity of this vertex and relationships between
the three solitons’ parameters can then be determined using a generalization
of Rankine-Hugoniot jump conditions for a family of conservation laws of the
modulation equations. In order to account for the multivalued nature of the Y
soliton while using (1+1)-dimensional modulation conservation laws, we
introduce _modified Rankine-Hugoniot jump conditions_ , where the densities
and fluxes of the two smaller legs of the Y soliton are summed. Notably, we
show that the family of modulation conservation laws for which the Y soliton
vertex satisfies these modified Rankine-Hugoniot jump conditions is _infinite_
, contrary to the typical setting of hyperbolic conservation laws.

Our results imply that the Y soliton can be understood as a modified _Whitham
shock_ [51], a recent discovery defined as a discontinuity in the modulation
equations that satisfies Rankine-Hugoniot jump conditions and corresponds to a
travelling wave solution to the full equation. Whereas in classical shock
theory, a shock is obtained by taking a zero viscosity limit of a travelling
wave solution to a dissipative equation, a Whitham shock is obtained by taking
the zero dispersion limit of a travelling wave in a dispersive equation.

The modulation theory approach features a number of additional novel elements.
A physically inspired way to study resonant soliton interactions is by
formulating non-soliton initial value problems. The descriptions of the full
dynamics of such problems have been limited to numerical simulations, because
analytical methods have focused on exact multi-soliton solutions to the KP
equation (1) [15, 34, 25, 52, 53, 26, 18]. The modulation approach has the
advantage of being able to describe the dynamical evolution of non-soliton
initial conditions, not only the long-time, locally steady state. In this
paper, we demonstrate this advantage by analytically solving for the evolution
of V-shaped initial conditions, i.e. the Mach reflection problem [38, 15, 25,
26, 53, 18]. In addition, to the best of our knowledge, this work is the first
fully two-dimensional application of Whitham modulation theory, which
typically has been restricted to one-dimensional equations or reductions.
Consequently, this method holds promise for studying soliton resonance or near
resonance in equations where multi-soliton solutions are not known.

The remainder of the paper is organized as follows. In section 2, we summarize
the exact solution forms of the KP line soliton and Y soliton solutions, and
we review their zero dispersion limits. Next, in section 3, the derivation of
an infinite number of soliton modulation equations in conservation form is
presented. In section 4, we show how the Y soliton satisfies (modified) jump
conditions for the infinite modulation conservation laws. These results are
applied to the classical Mach reflection problem in section 5, followed by a
discussion in section 6. Appendix A contains a proof relevant to our results.
In section 5, our analysis is quantitatively supported by a pseudospectral
numerical solver based on [17], the details of which are given in [49].

##  2 Exact soliton solutions to the KP equation

###  2.1 Exact soliton solutions

A large class of multi-soliton solutions of the KP equation (1) can be
expressed as [53, 54]:

$u(x,y,t)=12\,\frac{\partial^{2}}{\partial x^{2}}[\log\tau(x,y,t)]\,,$ (3a)
where the tau function $\tau(x,y,t)$, is given by the Wronskian of
exponentials
$\tau(x,y,t)=\mathop{\mathrm{Wr}}\nolimits(f_{1},\dots,f_{N})\,,\qquad
f_{n}(x,y,t)=\sum_{m=1}^{M}A_{n,m}e^{\theta_{m}(x,y,t)}\,,$ (3b) and the
phases $\theta_{1},\dots,\theta_{M}$ are
$\theta_{m}(x,y,t)=\frac{1}{\sqrt{12}\epsilon}\left(k_{m}x+\frac{k_{m}^{2}}{2}y-\frac{k_{m}^{3}}{3}t+\theta_{m,0}\right)\,.$
(3c)

The solution (3a) is uniquely determined by the phase parameters
$k_{1},\dots,k_{M}$ and the coefficient matrix $A=(A_{n,m})$, plus the
translation constants $\theta_{1,0},\dots,\theta_{M,0}$. Without loss of
generality, one can take the phase parameters to be strictly ordered so that
$k_{1}<\cdots<k_{M}$.

Generically, the above representation produces a solution with exactly $N$
asymptotic line solitons as $y\to\infty$ and $M-N$ asymptotic line solitons as
$y\to-\infty$ [54]. The amplitude and slope of each asymptotic line soliton
are completely determined by the phase parameters $k_{1},\dots,k_{M}$. In the
simplest nontrivial case, obtained when $N=1$ and $M=2$, one recovers the
soliton solution (2), shown in figure 2$(a)$. It is convenient to label the
two phase parameters as $k_{-}$, $k_{+}$, $k_{-}<k_{+}$ in this case. The
amplitude and slope parameters $a$ and $q$ are given by

$a=\frac{1}{4}(k_{+}-k_{-})^{2}\,,\qquad q=\frac{1}{2}(k_{+}+k_{-})\,.$ (4)

The inverse map to (4) is

$k_{-}=q-\sqrt{a}\,,\qquad k_{+}=q+\sqrt{a}\,.$ (5)

The next simplest case is obtained when $N=1$ and $M=3$ and represents the
resonant soliton solution [55], shown in figure 2$(b)$. As $y\to-\infty$ there
are two solitons with phase parameters $(k_{1},k_{2})$ and $(k_{2},k_{3})$,
respectively. Conversely, as $y\to\infty$ one observes only one soliton
specified by the phase parameters $(k_{1},k_{3})$. Consequently, the Y soliton
is determined (up to translation constants) by three phase parameters $k_{1}$,
$k_{2}$, and $k_{3}$. For leg $j$ of the Y soliton, $j\in\\{1,2,3\\}$,
counting clockwise from the bottom left (see figure 2$(b)$), the physical
variables $(a_{i},q_{i})$ are determined by the phase parameters as

$\displaystyle\begin{split}\begin{aligned}
&q_{1}=\frac{1}{2}(k_{2}+k_{1}),\qquad&&q_{2}=\frac{1}{2}(k_{3}+k_{1}),\qquad&&q_{3}=\frac{1}{2}(k_{3}+k_{2}),\\\
&a_{1}=\frac{1}{4}(k_{2}-k_{1})^{2},\qquad\qquad&&a_{2}=\frac{1}{4}(k_{3}-k_{1})^{2},\qquad\qquad&&a_{3}=\frac{1}{4}(k_{3}-k_{2})^{2}.\end{aligned}\end{split}$
(6)

Due to the ordering $k_{1}<k_{2}<k_{3}$, one can see that the top leg (leg 2,
also known as the Mach stem) is also the largest amplitude soliton. The
velocity $(s_{x},s_{y})$ of the Y soliton vertex is obtained by assuming a
traveling wave solution satisfying $\theta_{1}=\theta_{2}=\theta_{3}$ (3c),
which yields

$(s_{x},s_{y})=\left(-\frac{1}{3}(k_{1}k_{2}+k_{1}k_{3}+k_{2}k_{3}),\frac{2}{3}(k_{1}+k_{2}+k_{3})\right).$
(7)

Note that by applying the symmetry transformation $y\to-y$, one can apply the
same analysis to a Y soliton with two legs as $y\to\infty$ and one leg as
$y\to-\infty$. The expressions for the velocities (7) remain unchanged.
Although the Y soliton is a special exact solution to the KP equation (1),
Y-shaped interactions are ubiquitous as building blocks for higher-order
solutions to the KP equation [53].

Since the KP equation remains unchanged under the pseudorotation
transformation of the form [45, 56],

$q^{\prime}=Q+q,\qquad x^{\prime}=x+Qy-Q^{2}t,\qquad y^{\prime}=y-2Qt,\qquad
Q\in\mathds{R},$ (8)

one can determine $s_{x}$ in terms of the remaining Y soliton parameters. For
a general Y soliton with parameters (6), we consider the rotated image of the
Y soliton with parameters denoted by ′ and a vertical stem $q_{2}^{\prime}=0$.
The vertical velocity of the Y soliton image from (7) is
$s_{y}^{\prime}=2k_{2}^{\prime}/3$ and its horizontal velocity must be the
velocity of the stem $s_{x}^{\prime}=a_{2}/3$ (see (2)). If we now rotate the
image Y soliton to the initial Y soliton via a rotation of $Q=q_{2}$, we can
use the symmetries (8) to obtain the horizontal velocity that is equivalent to
(7),

$s_{x}=\frac{a_{2}}{3}+q_{2}^{2}-q_{2}s_{y}.$ (9)

###  2.2 Zero dispersion limit of KP solitons

We will study the Y soliton in the zero dispersion ($\epsilon\to 0$) limit.
The equation for the potential $\phi$ of a single line soliton (2) where
$u=\epsilon\phi_{x}$ is

$\phi(x,y,t)=\sqrt{12}q+\sqrt{12a}\,\tanh\left(\sqrt{\frac{a}{12}}\frac{x+qy-
ct}{\epsilon}\right),\qquad c=\frac{a}{3}+q^{2},$ (10)

where the first term is a constant of integration that is well-defined by
using the tau function form (3a). The zero dispersion parameter $\epsilon$
represents the transition width of the $\tanh$ function, which corresponds to
the width of the soliton (2). The limit of the potential (10) as $\epsilon\to
0$ becomes a discontinuous step,

$\lim_{\epsilon\to
0}\phi(x,y,t)=\sqrt{12}\left[(q-\sqrt{a})+2\sqrt{a}\,H\left(x+qy-
ct\right)\right],$ (11)

where $H$ is the Heaviside (step) function. Thus, the corresponding soliton
solution to the KP equation in the zero dispersion limit, which is the scaled
$x$-derivative of (11), is a zero-width impulse with parameters $a$ and $q$,
where $a$ is the impulse strength and $q$ is a measure of the slope of the
line in the $x$-$y$ plane. In other words, the soliton is localized along the
line $x+qy=ct$. This interpretation of a line soliton is similar to what has
been called a _soliton graph_ [53], with the zero dispersion limit
corresponding to the so-called “tropical” limit.

Consequently, in the zero dispersion limit, the three line solitons which
compose a Y soliton are localized to three semi-infinite lines with a single
point of intersection. This intersection point is a travelling discontinuity
between all three zero-width, constant line soliton impulses. It is the
behaviour of this point that we will interpret with Whitham modulation theory.

##  3 Soliton modulation theory

###  3.1 Derivation of one-dimensional soliton modulation equations

In this section, we will derive one-dimensional modulation conservation laws
for the parameters $(a,q)$ of the line soliton solution (2). These equations
describe the parameters’ evolution in the zero dispersion limit where the line
soliton reduces to the scaled derivative of (11), an impulse along a curve in
the $x$-$y$ plane. Fully two-dimensional soliton modulation equations have
been derived in prior work as [48, 47]

$\begin{bmatrix}a\\\
q\end{bmatrix}_{t}+\begin{bmatrix}\frac{1}{3}a-q^{2}&-\frac{4}{3}aq\\\
-\frac{1}{3}q&\frac{1}{3}a-q^{2}\end{bmatrix}\begin{bmatrix}a\\\
q\end{bmatrix}_{x}+\begin{bmatrix}2q&\frac{4}{3}a\\\
\frac{1}{3}&2q\end{bmatrix}\begin{bmatrix}a\\\ q\end{bmatrix}_{y}=0.$ (12)

We will use the following theorem, proved in Appendix A, to justify
considering a one-dimensional reduction of (12).

######  Theorem 3.1.

If $(a_{x},q_{x})$ and $(a_{y},q_{y})$ are well-defined and continuous, the
modulation equations (12) for the parameters $(a,q)$ of a KP line soliton (2)
can always be reduced to an equivalent (1+1)-dimensional system of equations.

Theorem 3.1 can be intuitively understood by the fact that since a line
soliton is localized to a curve in the zero dispersion limit, modulations in
one spatial dimension can be locally mapped to modulations in the other
spatial dimension using the line equation $x+qy=ct$. Consequently, in most of
this paper, without loss of generality we assume $x$-independence of the
modulation parameters. In section 4.4, we will show how to map these results
to $y$-independent parameters. As an aside, a natural approach for line
soliton modulations would be to use an arc length parametrization, similar to
equations obtained for modulations of line dark solitons of the nonlinear
Schrödinger equation [57].

In light of the above discussion, using the dispersion parameter $\epsilon$,
we assume the $x$-independent asymptotic expansion

$u(x,y,t;\epsilon)=u_{0}(\theta,y,t)+\epsilon u_{1}(\theta,y,t)+\cdots,$ (13)

with the leading order modulated soliton ansatz

$\begin{split}u_{0}(\theta,y,t;\epsilon)&=a(y,t)\operatorname{sech}^{2}\left(\sqrt{\frac{a(y,t)}{12}}\frac{\theta}{\epsilon}\right),\\\
\theta_{x}&=1,\quad\theta_{y}=q(y,t),\quad\theta_{t}=-\left(\frac{a(y,t)}{3}+q^{2}(y,t)\right),\end{split}$
(14)

and where every $u_{i}$ and their derivatives in (13) vanish as
$\theta\to\pm\infty$. Due to division by $\epsilon$ in (14), $\theta/\epsilon$
is a fast variable compared to $y$ and $t$. We require the compatibility
conditions $\theta_{xy}=\theta_{yx}$, $\theta_{xt}=\theta_{tx}$, and
$\theta_{yt}=\theta_{ty}$. The first two compatibility conditions are
satisfied by the assumption that $a$ and $q$ are $x$-independent, i.e.
$a_{x}\equiv q_{x}\equiv 0$. The third condition $\theta_{yt}=\theta_{ty}$
yields the equation

$\theta_{yt}-\theta_{ty}=q_{t}+\left(\frac{a}{3}+q^{2}\right)_{y}=0,$ (15)

which immediately gives us a modulation equation for $q$ in conservation form.
This equation (15) is known as conservation of waves.

We will now average the conservation laws of the KP equation over the fast
variable $\theta$ to obtain a number of additional modulation conservation
laws for the soliton amplitude [58, 40, 41]. Consider a general conservation
law

$\frac{\partial}{\partial t}F+\frac{\partial}{\partial
x}G+\frac{\partial}{\partial y}H=0,$ (16)

where $F,G,H$ depend on the derivatives of the potential $\phi$ such that
$u=\epsilon\phi_{x}$. Based on the modulated soliton form (14), we can expand
the derivatives as

$\frac{\partial}{\partial
t}\mapsto-\frac{c}{\epsilon}\frac{\partial}{\partial\theta}+\frac{\partial}{\partial
t},\qquad\frac{\partial}{\partial
x}\mapsto\frac{1}{\epsilon}\frac{\partial}{\partial\theta},\qquad\frac{\partial}{\partial
y}\mapsto\frac{q}{\epsilon}\frac{\partial}{\partial\theta}+\frac{\partial}{\partial
y}.$ (17)

After inserting the derivatives (17) and the expansion (13) into (16), we then
apply the integral operator
$\bar{F}=\int_{-\infty}^{\infty}F(u)\,\mathrm{d}\theta$ to each element in the
equation. All total derivatives with respect to $\theta$ disappear, since we
assumed that every $u_{i}$ and their derivatives vanish as
$\theta\to\pm\infty$. We are therefore left with

$\int_{-\infty}^{\infty}\left(\frac{\partial F}{\partial t}+\frac{\partial
G}{\partial x}+\frac{\partial H}{\partial
y}\right)\mathrm{d}\theta=\frac{\partial\overline{F}}{\partial
t}+\frac{\partial\overline{H}}{\partial y}+\mathcal{O}(\epsilon).$ (18)

As a consequence of integrability, the KP equation (1) admits an infinite
number of conservation laws of the form (16). We list three of them here [59,
60], given in terms of the potential $\phi$ where $u=\epsilon\phi_{x}$:

$\displaystyle\left(\phi_{x}^{2}\right)_{t}+\left(2\epsilon^{2}(\phi_{x}\phi_{xxx}-\phi_{xx}^{2})+\frac{2\epsilon}{3}\phi_{x}^{3}-\phi_{y}^{2}\right)_{x}+\left(2\phi_{x}\phi_{y}\right)_{y}$
$\displaystyle=0,$ (19a)
$\displaystyle\begin{split}\left(\phi_{x}\phi_{y}\right)_{t}+\left(2\epsilon^{2}(\phi_{y}\phi_{xxx}-\phi_{xx}\phi_{xy})+\epsilon\phi_{y}\phi_{x}^{2}+\phi_{t}\phi_{y}\right)_{x}+\\\
+\left(\epsilon^{2}\phi_{xx}^{2}+\phi_{y}^{2}-\frac{\epsilon}{3}\phi_{x}^{3}-\phi_{t}\phi_{x}\right)_{y}&=0,\end{split}$
(19b)
$\displaystyle\left(\frac{\epsilon\phi_{x}^{3}}{3}-\epsilon^{2}\phi_{xx}^{2}+\phi_{y}^{2}\right)_{t}+\left(2\epsilon^{2}(\phi_{xt}\phi_{xx}-\phi_{t}\phi_{xxx})-\epsilon\phi_{t}\phi_{x}^{2}-\phi_{t}^{2}\right)_{x}-\left(2\phi_{t}\phi_{y}\right)_{y}$
$\displaystyle=0.$ (19c)

These equations correspond to conservation of momentum in $x$ (19a), momentum
in $y$ (19b), and energy (19c) in shallow water applications. After
transforming the derivatives according to (17), inserting the expansion (13),
and applying the integral operator, the zero dispersion limit of (19) becomes

$\displaystyle\left(\overline{u_{0}^{2}}\right)_{t}+\left(2q\overline{u_{0}^{2}}\right)_{y}$
$\displaystyle=0$ (20a)
$\displaystyle\left(q\overline{u_{0}^{2}}\right)_{t}+\left(\overline{u_{0\theta}^{2}}-\frac{\overline{u_{0}^{3}}}{3}+(q^{2}+c)\overline{u_{0}^{2}}\right)_{y}$
$\displaystyle=0$ (20b)
$\displaystyle\left(-\overline{u_{0\theta}^{2}}+\frac{\overline{u_{0}^{3}}}{3}+q^{2}\overline{u_{0}^{2}}\right)_{t}+\left(2qc\overline{u_{0}^{2}}\right)_{y}$
$\displaystyle=0.$ (20c)

Upon substituting (14), the integrals in (20) can be calculated explicitly to
yield the modulation conservation laws

$\displaystyle\left(a^{3/2}\right)_{t}+\left(2qa^{3/2}\right)_{y}$
$\displaystyle=0,$ (21a)
$\displaystyle\left(qa^{3/2}\right)_{t}+\left(\frac{2}{15}a^{5/2}+2q^{2}a^{3/2}\right)_{y}$
$\displaystyle=0,$ (21b)
$\displaystyle\left(\frac{a^{5/2}}{5}+q^{2}a^{3/2}\right)_{t}+\left(\frac{2}{3}qa^{5/2}+2q^{3}a^{3/2}\right)_{y}$
$\displaystyle=0.$ (21c)

The expression (21a) corresponds to conservation of wave action for soliton
modulations. Expanding the derivatives in each of (21) yields equivalent
modulation equations for $a$ and $q$. Together with (15), the $x$-independent
modulation equations for $a$ and $q$ can then be written in matrix form as the
$x$-independent reduction of (12),

$\begin{bmatrix}a\\\ q\end{bmatrix}_{t}+\begin{bmatrix}2q&\frac{4}{3}a\\\
\frac{1}{3}&2q\end{bmatrix}\begin{bmatrix}a\\\ q\end{bmatrix}_{y}=0.$ (22)

Importantly, all three modulation conservation laws (21) are consistent with
the system (22). This system (22) has previously been derived using multiple
scales [61, 49] and averaged Lagrangian methods [47]. The KP equation has an
additional conservation law, conservation of mass (with density
$u=\epsilon\phi_{x}$), but we do not consider this law as it yields results
inconsistent with (22) and (21). Note also that the modulation system (22) and
therefore all modulation conservation laws [62] possess the scaling symmetry

$(y,t,a,q)\longrightarrow(\lambda^{2}y^{\prime},\lambda^{3}t^{\prime},\lambda^{2}a^{\prime},\lambda
q^{\prime}).$ (23)

The infinite number of KP conservation laws suggests that, in principle, one
could continue this averaging process to find an infinite number of modulation
conservation laws. In the next section we will use a more direct approach to
obtain an infinite number of modulation conservation laws.

###  3.2 Derivation of generalized conservation laws

In this section we generate an infinite family of conservation laws consistent
with the modulation system (22), following an approach developed by G.B.
Whitham [40]. A general soliton modulation conservation law in $y$, like those
given in (21), has the form

$\frac{\partial f}{\partial t}+\frac{\partial h}{\partial y}=0,$ (24)

where $f=f(a,q)$ and $h=h(a,q)$. Applying the chain rule and inserting the
evolution equations (22) yields

$\left(-2qf_{a}-\frac{1}{3}f_{q}+h_{a}\right)a_{y}+\left(-\frac{4}{3}af_{a}-2qf_{q}+h_{q}\right)q_{y}=0.$
(25)

Since (25) must be true for all values of $q_{y}$ and $a_{y}$, the terms in
the parentheses must both equal zero, leading to two compatibility equations
for $f$ and $h$. Let us consider a general density of the form

$f(a,q)=\sum^{n}_{i=1}c_{i}a^{i+\frac{1}{2}}q^{2(n-i)},\qquad c_{1}=1.$ (26)

This density ansatz is motivated by analogy with (21a) (which corresponds to
(26) with $n=1$) and (21c) (which corresponds to (26) with $n=2$). In
addition, the ansatz (26) satisfies the scaling symmetries (23) of the
modulation system (22). Taking the partial derivatives of (26) and inserting
into the compatibility conditions from (25) yields

$\displaystyle h_{a}$
$\displaystyle=\sum^{n}_{i=1}\left[(2i+1)q+\frac{2}{3}(n-i)\right]c_{i}a^{i+\frac{1}{2}}q^{2(n-i+\frac{1}{2})},$
(27a) $\displaystyle h_{q}$
$\displaystyle=\sum^{n}_{i=1}\left[\frac{2}{3}(2i+1)+4(n-i)\right]c_{i}a^{i+\frac{1}{2}}q^{2(n-i)}.$
(27b)

Requiring equality of coefficients for $h_{aq}=h_{qa}$ gives a recursion
relation for $c_{i}$ in terms of $c_{i-1}$, which can be rewritten explicitly
as:

$c_{i}=\frac{3i}{(n-i)(2n)(2n-1)}{2n\choose 2i+1},\quad i=1,\dots,n-1,\qquad
c_{n}=\frac{3}{(2n+1)(2n-1)}.$ (28)

By integrating $h_{q}$ (27b) with respect to $q$ we obtain the flux as

$h(a,q)=\sum^{n}_{i=1}d_{i}a^{i+\frac{1}{2}}q^{2(n-i)+1},$ (29)

where the coefficients $d_{i}$ are defined as

$d_{i}=\left(2-\frac{4(i-1)}{3(2i-2n-1)}\right)c_{i},\qquad i=1,\dots,n.$ (30)

Up to an arbitrary additive constant, the same result (29) is also obtained by
integrating (27a). Since the constant is arbitrary, we set it to zero. For
$n=1,2,\dots$, the equations (24), (26), (28), (29), and (30) together form an
infinite family of modulation conservation laws that are consistent with the
modulation system (22). This family includes (21a) and (21c) as the $n=1$ and
$n=2$ cases, respectively. The $n=3$ case is

$\left(\frac{3}{35}a^{7/2}+\frac{6}{5}q^{2}a^{5/2}+q^{4}a^{3/2}\right)_{t}+\left(\frac{2}{5}qa^{7/2}+\frac{44}{15}q^{3}a^{5/2}+2q^{5}a^{3/2}\right)_{y}=0.$
(31)

We can use a different density ansatz to obtain a second infinite family of
modulation conservation laws. This family includes (21b) as the $n=1$ case.
The process is the same as above, so we simply state the results. The density
$\tilde{f}(a,q)$ takes the form

$\displaystyle\tilde{f}(a,q)$
$\displaystyle=\sum^{n}_{i=1}\tilde{c}_{i}a^{i+\frac{1}{2}}q^{2(n-i)+1},\qquad\quad\tilde{c}_{1}=1,$
(32a) $\displaystyle\tilde{c}_{i}$
$\displaystyle=\frac{(2i-2n-3)(i-n-1)}{(i-1)(2i+1)}\tilde{c}_{i-1},\qquad
i=2,\dots,n.$ (32b)

The corresponding flux $\tilde{h}(a,q)$ can be found to be

$\displaystyle\tilde{h}(a,q)$
$\displaystyle=\sum^{n+1}_{i=1}\tilde{d}_{i}a^{i+\frac{1}{2}}q^{2(n-i)+2},\qquad\tilde{d}_{1}=2,\qquad\tilde{d}_{n+1}=\frac{2}{9+6n}\tilde{c}_{n},$
(33a) $\displaystyle\tilde{d}_{i}$
$\displaystyle=2\tilde{c}_{i}+\frac{2n-2i+1}{3(i+3/2)}\tilde{c}_{i-1},\qquad
i=2,\dots,n.$ (33b)

##  4 Shock solutions of the soliton modulation equations

In this section we show that the Y soliton can be understood as a
discontinuous, weak solution to all the modulation conservation laws given in
section 3. In particular, the velocity of the vertex $s_{y}$ of the Y soliton
and each legs’ parameters satisfy Rankine-Hugoniot conditions for the
conservation of waves equation (15) and _modified_ Rankine-Hugoniot jump
conditions for all other modulation conservation laws given above. We also
summarize some properties of the Y soliton shock.

###  4.1 Shock solution to general KP conservation laws

We can use jump conditions arising from any averaged KP conservation law to
determine the correct velocity of the Y soliton. However, it is evident from
figure 2$(b)$ that the modulation equations for $a(y,t)$ and $q(y,t)$ are
_multivalued_ in $y$; for example, for $y<z(t)$, $\dot{z}=s_{y}$, $a(y,t)$
takes on two different constant values: $a(y,t)=a_{1}$ and $a(y,t)=a_{3}$.
Consequently, the calculation of the jump condition must be modified from
traditional Rankine-Hugoniot jump conditions, which are defined using the
difference operator $\llbracket u\rrbracket=u_{+}-u_{-}$ between two states
$u_{\pm}$. Instead, we define a modified difference operator:

######  Definition 4.1.

Consider a multivalued function $u(y)$ with a discontinuity at $y=z(t)$ and

$\lim_{y\to z(t)^{+}}u=\\{u_{+}^{(1)},\dots,u_{+}^{(n)}\\},\qquad\lim_{y\to
z(t)^{-}}u=\\{u_{-}^{(1)},\dots,u_{-}^{(m)}\\},\qquad\dot{z}(t)=s_{y}(t).$

The _multivalued difference operator_ $\llbracket\,\rrbracket_{n,m}$ is
defined as

$\llbracket
u\rrbracket_{n,m}=\sum_{j=1}^{n}u_{+}^{(j)}-\sum_{j=1}^{m}u_{-}^{(j)}.$

Using the multivalued difference operator, we can define modified Rankine-
Hugoniot conditions:

######  Definition 4.2.

For a conservation law $(f(u))_{t}+(h(u))_{y}=0$, where $u(y,t)$ is
multivalued and discontinuous as in definition 4.1, the discontinuity with
velocity $s_{y}$ is calculated using _modified Rankine-Hugoniot conditions_ as

$s_{y}\llbracket f(u)\rrbracket_{n,m}=\llbracket h(u)\rrbracket_{n,m}\,.$

Since both $a(y,t)$ and $q(y,t)$ are multivalued for either $y<z(t)$ or
$y>z(t)$ for the Y soliton, we will apply modified Rankine-Hugoniot jump
conditions to the modulation conservation laws derived above. The summation of
multivalued densities and fluxes in definition 4.2 will be justified below.
First, we demonstrate how to perform this calculation using the wave action
conservation law (21a), for which modified Rankine-Hugoniot jump conditions
take the form

$s_{y}=\frac{\left\llbracket 2qa^{3/2}\right\rrbracket_{n,m}}{\left\llbracket
a^{3/2}\right\rrbracket_{n,m}}\,.$ (34)

The multivalued difference operator $\llbracket\,\rrbracket_{n,m}$ is applied
by adding the appropriate quantities that share the same half plane in $y$, in
this case legs 1 and 3 (see figure 2$(b)$), so that $n=1$, $m=2$, and

$s_{y}=\frac{\left(2q_{2}a_{2}^{3/2}\right)-\left(2q_{1}a_{1}^{3/2}+2q_{3}a_{3}^{3/2}\right)}{\left(a_{2}^{3/2}\right)-\left(a_{1}^{3/2}+a_{3}^{3/2}\right)}=\frac{2}{3}\left(k_{1}+k_{2}+k_{3}\right),$
(35)

where the simplification from physical to phase variables follows from (6)
after some algebra. The result (35) agrees with the known velocity $s_{y}$ of
the Y soliton vertex in $y$ (7). We remark that modified jump conditions for
the zero dispersion limit of the Y soliton applied to any of the modulation
conservation equations (21) will also correctly determine $s_{y}$.

One can justify the adding of the fluxes and densities in $y$ as follows.
Consider any conservation law of the full KP equation of the form (16). A
general travelling wave solution that is localized in $x$ with velocity
$(s_{x},s_{y})$ has the form

$u=u(\xi,\eta),\qquad\xi=x-s_{x}t,\qquad\eta=y-s_{y}t.$ (36)

Inserting (36) into (16) yields

$-s_{x}F_{\xi}-s_{y}F_{\eta}+G_{\xi}+H_{\eta}=0.$ (37)

We now fix $\eta$ and integrate (37) with respect to $\xi$ from $-\infty$ to
$\infty$. This definite integral can be calculated for the first and third
terms of (37) by evaluating $F$ and $G$ at $\xi\to\pm\infty$, where they are
zero due to the localized nature of $u$. Thus, we obtain

$\left(\int_{-\infty}^{\infty}-s_{y}F(\xi,\eta)+H(\xi,\eta)\,\mathrm{d}\xi\right)_{\eta}=0.$

This expression can be integrated with respect to $\eta$ to yield

$-s_{y}\overline{F}(\eta)+\overline{H}(\eta)=A,$ (38)

where $A$ is a constant, and the average, denoted e.g. by $\overline{F}$, is
calculated by integrating over the real line in the $\xi$-direction.
Evaluating (38) at two different values $\eta^{+}\gg 1$ and $-\eta^{-}\gg 1$
allows us to eliminate the constant $A$, and we can solve for $s_{y}$ to
obtain

$s_{y}=\frac{\overline{H}(\eta^{+})-\overline{H}(\eta^{-})}{\overline{F}(\eta^{+})-\overline{F}(\eta^{-})},$

which has an identical form to (34). For the Y soliton given in section 2,
averaging the flux and density over $\xi$ from $-\infty$ to $\infty$, at some
$-\eta^{-}\gg 1$, is equivalent to adding together the fluxes and densities of
the two smaller legs (legs 1 and 3). Thus, $\overline{F}(\eta^{-})$ can be
expressed as $\overline{F}_{1}+\overline{F}_{3}$, where $\overline{F}_{i}$
denotes the value of the averaged flux $\overline{F}$ at the $i^{\rm th}$ leg.
Similarly, $\overline{H}(\eta^{-})=\overline{H}_{1}+\overline{H}_{3}$. Our
final expression for the velocity in $y$ becomes

$s_{y}=\frac{\overline{H}_{2}-\left(\,\overline{H}_{1}+\overline{H}_{3}\right)}{\overline{F}_{2}-\left(\,\overline{F}_{1}+\overline{F}_{3}\right)},$

which corresponds to the addition of fluxes and densities in the averaged
conservation law in (35). In the zero dispersion limit, the Y soliton vertex
region shrinks to a point, and we merely require $\eta^{+}>0$ and
$\eta^{-}<0$. A similar construction holds if the multivalued pair of Y
soliton legs occurs for $y>z(t)$ where now $n=2$, $m=1$ in (34).

###  4.2 Shock solutions to families of conservation laws

In this section, we will show that the Y soliton satisfies modified Rankine-
Hugoniot conditions (see definition 4.2) for the infinite family of
conservation laws given above. Specifically, we will prove the following:

######  Proposition 4.1.

Consider a modulation conservation law (24) with density given by (26) and
(28) and flux given by (29) and (30), respectively. Then modified Rankine-
Hugoniot jump conditions for the conservation law applied to the Y soliton
parameters (6) yield

$\frac{\left\llbracket h(a,q)\right\rrbracket_{1,2}}{\left\llbracket
f(a,q)\right\rrbracket_{1,2}}=\frac{2}{3}(k_{1}+k_{2}+k_{3})=s_{y}.$ (39)

###### Proof.

We proceed by rewriting the density for a single soliton (26) in terms of
phase variables $s=q+\sqrt{a}$ and $r=q-\sqrt{a}$ (5) as

$f(s,r)=\sum_{i=1}^{n}\frac{c_{i}}{2^{2n+1}}(s-r)^{2i+1}(s+r)^{2(n-i)}.$ (40)

The density polynomial (40) can be rewritten as

$f(s,r)=(s^{2n+1}-r^{2n+1})\sum_{i=1}^{n}\frac{c_{i}}{2^{2n+1}}+(r^{2n}s-rs^{2n})\sum_{i=1}^{n}\frac{c_{i}}{2^{2n+1}}(4i-2n+1).$

These terms come from expanding the binomials in (40) and distributing for the
product. Any term that contains a power of less than $2n$ cancels. Similarly,
the flux (29) can be rewritten as

$\displaystyle h(s,r)$
$\displaystyle=\sum_{i=1}^{n}\frac{d_{i}}{2^{2n+1}}(s-r)^{2i+2}(s+r)^{2(n-i)+1},$
$\displaystyle\begin{split}&=(s^{2n+2}-r^{2n+2})\sum_{i=1}^{n}\frac{d_{i}}{2^{2n+2}}+(r^{2n}s-rs^{2n})\sum_{i=1}^{n}\frac{d_{i}}{2^{2n+2}}(4i-2n)+\\\
&\quad+(r^{2n-1}s^{2}-r^{2}s^{2n-1})\sum_{i=1}^{n}\frac{d_{i}}{2^{2n+2}}(2n^{2}-n-8in+8i^{2}-1).\end{split}$

If we now insert the phase variables $k_{i}$ for each leg, evaluate
$\left\llbracket f\right\rrbracket_{1,2}$, and multiply by the velocity
$s_{y}=\frac{2}{3}(k_{1}+k_{2}+k_{3})$, we find after some algebra that this
product is equal to $\left\llbracket h\right\rrbracket_{1,2}$ if and only if
the following equalities hold

$\sum_{i=1}^{n}\frac{c_{i}}{3}(4i-2n+1)=\sum_{i=1}^{n}\frac{d_{i}}{2}(2i-n)=\sum_{i=1}^{n}-\frac{d_{i}}{4}(2n^{2}-n-8in+8i^{2}-1).$
(41)

Thus, if we can prove (41), we have shown (39). Using (28) and (30), a
calculation demonstrates that (41) holds. ∎

Thus, the parameters of the resonant Y soliton satisfy modified Rankine-
Hugoniot jump conditions for an infinite family of conservation laws in $y$.
We omit the similar proof for the second family of conservation laws given in
(32) and (33).

###  4.3 Shock solution to conservation of waves

In this section, we use _ordinary_ Rankine-Hugoniot jump conditions for the
modulation equation for $q$ (15) that arises from conservation of waves.
Instead of modified Rankine-Hugoniot jump conditions, we choose any two legs
of the Y soliton and insert them into the ordinary jump conditions for (15).
The justification for using modified Rankine-Hugoniot jump conditions above
depended on the fact that the modulation conservation laws were derived from
KP conservation laws. In contrast, the wave action conservation law (15) is a
kinematic condition arising from the definition of our modulation ansatz (14).
A similar kinematic calculation to determine the speed of multisoliton
structures was utilized in [63].

Without loss of generality, we consider the jump between legs 1 and 2 (see
figure 2$(b)$). Then, ordinary jump conditions
($\llbracket\,\rrbracket\equiv\llbracket\,\rrbracket_{1,1}$) for (15) become
(with simplification using (6))

$\displaystyle s_{y}\llbracket q\rrbracket$
$\displaystyle=\left\llbracket\frac{a}{3}+q^{2}\right\rrbracket,$ (42)
$\displaystyle s_{y}$
$\displaystyle=\frac{\frac{1}{3}(k_{2}-k_{1})^{2}+(k_{2}+k_{1})^{2}-\frac{1}{3}(k_{3}-k_{1})^{2}-(k_{3}+k_{1})^{2}}{2((k_{2}+k_{1})-(k_{3}+k_{1}))}=\frac{2}{3}(k_{1}+k_{2}+k_{3}),$
(43)

which is the correct velocity $s_{y}$ (7). One can verify that the jump
condition (42) yields the correct velocities for the choice of _any two legs_
from the three that make up the resonant Y soliton. This is true even for the
two legs in the multivalued region, i.e. inserting legs 1 and 3 into the jump
conditions in $y$ (42).

We justify this result as follows. The velocity of the vertex of the Y soliton
is calculated directly by setting $\theta_{1}=\theta_{2}=\theta_{3}$, where
$\theta_{i}$ is the non-modulated phase (2) of the $i^{\rm th}$ soliton leg.
Consider $\theta_{i}=\theta_{j}$, which can be rewritten as

$(q_{i}-q_{j})y=\left(\frac{a_{i}}{3}+q_{i}^{2}-\frac{a_{j}}{3}+q_{j}^{2}\right)t,$

up to an overall additive phase constant that we have set to zero. Since
$y/t=s_{y}$, the above equation is precisely the jump condition (42). Thus,
the requirement that $\theta_{i}=\theta_{j}$ is equivalent to calculating the
shock velocity $s_{y}$ using ordinary jump conditions. The jump conditions
(42) are essentially kinematic conditions ensuring that the three solitons are
travelling together.

###  4.4 Shocks in $y$-independent variables

In this section, we briefly show how the results from the above sections using
$x$-independent equations can be repeated for $y$-independent equations,
correctly obtaining the velocity $s_{x}$ of the Y soliton vertex in $x$. In
order to transform $x$-independent modulation equations (22) to
$y$-independent modulation equations, one can use the identification

$\partial_{y}\to-q\partial_{x},\qquad\partial_{t}\to\partial_{t}+c\partial_{x},$
(44)

where the first expression in (44) transforms the derivative with respect to
$y$ to a derivative with respect to $x$ (assuming $q\neq 0$). The second
expression in (44) reflects that modulations in $x$ must occur at the velocity
$c$ of the soliton in $x$.

Assuming a modulation conservation law in $y$ with the general form
$f_{t}+h_{y}=0$, one can utilize the transformation (44) combined with
conservation of waves $q_{t}+c_{y}=0$ (c.f. (15)) to verify that there is an
equivalent modulation conservation law in $x$ having the form

$\left(\frac{f}{q}\right)_{t}+\left(\frac{fc}{q}-h\right)_{x}=0.$ (45)

With (45), one can extend all of the above results in $y$ to results in $x$.

###  4.5 Properties of the Y soliton shock solution

####  4.5.1 Construction and uniqueness of Y soliton shocks

In this section we will discuss how to construct the Y soliton using
modulation conservation laws and make some remarks regarding uniqueness. For
three line solitons travelling together there are eight free variables: six
soliton parameters and the velocity $(s_{x},s_{y})$. As the Y soliton is a
three-parameter family of solutions to the KP equation (1), at least five
equations are required to completely determine the Y soliton. The speed
relation (9) and two jump conditions for conservation of waves (15) yield
three independent equations. We obtain additional equations using the modified
jump conditions in definition 4.2 applied to (21) and the infinite family of
conservation laws. The eight unknown parameters of the Y soliton solve this
infinite system of algebraic equations.

A discontinuous solution to the modulation conservation laws as described
above is said to be _admissible_ if it corresponds to the zero dispersion
limit of a travelling wave solution to the KP equation (1). The Y soliton is
the only KP travelling wave solution consisting of three solitons joined
together, so any other discontinuous solution to the modulation conservation
laws with three solitons is therefore inadmissible. It remains an open
question as to how many modulation conservation laws are necessary to uniquely
determine the Y soliton solution, i.e. what is the minimal subset of the
infinite set of algebraic equations that does not yield inadmissible
solutions.

The Hilbert basis theorem implies that a finite subset of equations exists
that gives necessary and sufficient conditions to determine the unique
solution to the entire infinite set of equations (see e.g. [64]). The
simplest, smallest candidate for this subset consists of (9), two equations
from jump conditions for conservation of waves (15), and modified jump
conditions for conservation of momentum in $x$ (21a) and $y$ (21b). Numerical
calculations reveal that these five equations have a second (irrational)
solution family that is inadmissible. Including the modified jump condition
for conservation of energy (21c) does not eliminate this inadmissible family,
but the next conservation law (31) does eliminate it. Whether or not this set
of six algebraic equations is sufficient to guarantee a unique solution
corresponding to the admissible Y soliton is a question left for further
study.

Using a similar approach, we can show that _no two half solitons_ can be
connected by a single shock satisfying modified Rankine-Hugoniot jump
conditions. First, no KP travelling wave solution consisting of two half
solitons joined together exists, other than a soliton (trivially), so no
corresponding discontinuous solution can be admissible. Second, attempting to
solve for the parameters of such a solution using conservation laws yields a
combination of parameters with nonzero imaginary parts or negative amplitude.
Consider a shock solution to the modulation system (22) consisting of two half
solitons

$a(y,t)=\begin{cases}a_{2}&y>s_{y}t\\\ a_{1}&y<s_{y}t\end{cases},\qquad
q(y,t)=\begin{cases}q_{2}&y>s_{y}t\\\ q_{1}&y<s_{y}t\end{cases},$ (46)

where $a_{1},a_{2}>0$. In order for a shock to exist in the modulation
variables, ordinary jump conditions arising from conservation of waves (42),
as well as modified jump conditions for the conservation laws (21a) and (21b)
should all equal $s_{y}$. Setting these three jump conditions equal to each
other eliminates $s_{y}$ and gives two polynomial equations in $q$ and
$\sqrt{a}$. Given some lower soliton with parameters $(a_{1},q_{1})$, we
obtain a nonlinear system with two equations and two unknowns. Solving for
$(a_{2},q_{2})$ yields the following solutions for $a_{2}$

$a_{2}\in\bigg{\\{}\pm
a_{1},\left(-\frac{1}{2}\pm\frac{i}{2}\pm\sqrt{-1\pm\frac{i}{2}}\right)a_{1}\bigg{\\}},$
(47)

where each $\pm$ should be understood to operate independently from the
others, yielding a total of ten unique solutions for $a_{2}$. Given that
$a_{1}>0$, each solution for $a_{2}$ in (47) is either negative or complex and
is therefore not admissible. Consequently, it is impossible for initial
conditions as given in (46) to satisfy the modulation jump conditions unless
$a_{1}=a_{2}$, $q_{1}=q_{2}$. Any shock solution to the modulation system (22)
must therefore consist of at least three half solitons joined together.

####  4.5.2 Y soliton as shock in Riemann invariants

We find some additional properties of the Y soliton modulation shock by
examining the $x$-independent modulation system (22), which can be written in
diagonal form as [49]

$\frac{\partial}{\partial t}R_{\pm}+V_{\pm}\frac{\partial}{\partial
y}R_{\pm}=0,\qquad R_{\pm}=q\pm\sqrt{a},\qquad
V_{\pm}=2q\pm\frac{2}{3}\sqrt{a}=\frac{4}{3}R_{\pm}+\frac{2}{3}R_{\mp}\,,$
(48)

where $R_{\pm}$ are Riemann invariants and $V_{\pm}$ are the corresponding
characteristic velocities. As pointed out in [50], $R_{\pm}$ above are
identical to the phase variables $k\pm$ for a soliton (5). As a result, for
the Y soliton, one Riemann invariant is conserved for the jump from the Y
soliton stem (leg 2 in figure 2$(b)$) to the other legs (leg 1 or 3). From leg
2 to leg 3, $R_{+}=k_{3}$ is conserved. Examining characteristic velocities
(48) shows that this is a _compressive_ shock in the slow characteristic
family (i.e. $V_{-}(y>s_{y}t)<V_{-}(y<s_{y}t)$). From leg 2 to leg 1 in figure
2$(b)$, $R_{-}=k_{1}$ is conserved, which in this case is an _expansive_ shock
in the fast characteristic family (i.e. $V_{+}(y>s_{y}t)>V_{+}(y<s_{y}t)$).
For the jump from leg 1 to leg 3, neither Riemann invariant is conserved. This
jump in $x$ passes through $q=0$, where the $y$-independent modulation system
is nonstrictly hyperbolic, and in fact the $y$-independent soliton ansatz
cannot be defined for $q=0$ (cf. [50]).

The conservation of Riemann invariants gives us one more way to calculate the
velocity of the Y soliton $s_{y}$. For the jump where $R_{-}=k_{1}$ is
conserved, from leg 1 to leg 2, we can use ordinary jump conditions for the
hyperbolic equation for the Riemann invariant in $y$ (48) to correctly
calculate $s_{y}$, ignoring leg 3. A similar calculation holds for
conservation of $R_{+}$ from leg 2 to leg 3.

##  5 Application: V-shape initial conditions

The previous sections showed that the Y soliton can be understood as a shock
solution to an infinite family of modulation conservation laws. In this
section we apply these insights to V-shape initial conditions, i.e. the Mach
reflection problem, which in modulation variables is written as

$a(y,0)=a_{0},\qquad q(y,0)=\begin{cases}-q_{0}&y>0\\\
q_{0}&y<0\end{cases},\quad$ (49)

where $q_{0}>0$. The initial conditions (49) model the symmetrical oblique
interaction of two line solitons, which is identical to a single soliton
encountering an inward oblique corner through a method of images argument [25,
26] (see schematic in figure 1). This problem is well-studied [15, 25, 26, 53,
18], but the following modulation theory approach offers new insight on the
dynamics of the system evolution. Given a modulation solution $a(y,t)$ and
$q(y,t)$ to (49), the solution of the full KP equation (1) can be
reconstructed via

$u(x,y,t)=a(y,t)\,\mathrm{sech}^{2}\left(\sqrt{\frac{a(y,t)}{12}}\theta\right),\qquad\theta=x+\int_{0}^{y}q(y^{\prime},t)\,\mathrm{d}y^{\prime}-\int_{0}^{t}c(0,t^{\prime})\,\mathrm{d}t^{\prime},$
(50)

where the soliton velocity satisfies $c(y,t)=a(y,t)/3+q(y,t)^{2}$ (cf. (2),
(14)). Numerical evolutions of the full KP equation (1) with $\epsilon=1$ and
initial conditions (49) projected onto (50) (with smoothing of $q(y,0)$ to
prevent Gibbs phenomenon) are shown in figure 3 for $\sqrt{a_{0}}>q_{0}$ and
in figure 5 for $\sqrt{a_{0}}<q_{0}$. Note that the scale-free initial
condition (49) and modulation equations (12) imply that the $\epsilon\to 0$
limit for $t=\mathcal{O}(1)$ is equivalent to the $t\to\infty$ limit with
$\epsilon=\mathcal{O}(1)$.

Figure 3: Numerical evolution of V-shape initial conditions
$\eqref{eq:ArrowheadIC}$ projected onto (50) for the KP equation (1) with
$\epsilon=1$ for $a_{0}=1$ and $q_{0}=.7$ on $t\in\\{0,100,250\\}$.

###  5.1 Mach reflection

By directly examining the characteristic velocities (48), it is evident that
the solution to (49) must give rise to shocks in the modulation variables. For
$y>0$ in (49) the characteristic velocities are $V_{\pm}(y>0)=-2q_{0}\pm
2\sqrt{a_{0}}/3$, while for $y<0$, $V_{\pm}(y<0)=2q_{0}\pm 2\sqrt{a_{0}}/3$.
Since $q_{0}>0$, $V_{\pm}(y<0)>V_{\pm}(y>0)$, so both characteristic
velocities are compressive and should give rise to shocks in the modulation
equations. As Y solitons are shocks in modulation variables, we will look for
solutions to (49) consisting of two Y solitons.

Initial data (49) for the $x$-independent modulation equations (15) and (20)
constitute a Riemann problem for a hyperbolic system of two equations.
Traditionally, such problems are solved in terms of two waves—one slow and the
other fast—that are separated by a constant region. Each wave is either a
shock or rarefaction wave. Here, the multivalued structure of the Y shock and
the modified Rankine-Hugoniot conditions lead us to solve this Riemann problem
with four waves separated by constant regions. The four waves are visible in
the KP simulation of figure 3. We observe the two reflected Y shocks separated
by a constant stem $(a_{\rm i},q_{\rm i})$ and two rarefaction waves that lead
the growth of two new partial soliton legs from $a=0$ to $(a,q)=(a_{\rm
n},q_{\rm n})$ (upper left leg) and $(a,q)=(a_{\rm n},-q_{\rm n})$ (lower left
leg). The new legs arise as a consequence of the fact that modulation shocks
must be composed of at least three soliton legs to satisfy the jump conditions
(see section 4.5).

Figure 4: Comparison between analytical solution (dashed) with numerical
solution (solid) for v-shape initial conditions from figure 3 for $a$ (left)
and $q$ (right) on $t\in\\{0,50,150,250\\}$. The top panel shows the evolution
of the upper left leg as modelled by (53), and the bottom panel shows the
evolution of the original legs and the shock stem as described by (51). In all
cases the evolution approaches the long-time prediction (dashed-dotted). The
top plot includes a phase shift of $y_{0}=12$ in order to account for higher-
order effects.

We can determine the parameters of the stem $(a_{\rm i},q_{\rm i})$ by looking
for a state that shares one phase variable (Riemann invariant) with the top
and the bottom, solving the equations

$-q_{0}-\sqrt{a_{0}}=q_{\rm i}-\sqrt{a_{\rm i}},\qquad
q_{0}+\sqrt{a_{0}}=q_{\rm i}+\sqrt{a_{\rm i}}.$

We obtain that $q_{\rm i}=0$ and $\sqrt{a_{\rm i}}=\sqrt{a_{0}}+q_{0}$. From
the requirement that the physical variables of a Y soliton conserve phase
parameters (see section 4.5), one can determine that the new upper left leg
has parameters $(a_{\rm n},q_{\rm n})=(q_{0}^{2},\sqrt{a_{0}})$, with
symmetrical results for the bottom new leg. Using the jump conditions
calculated throughout this paper, we also determine that the velocity of the
top discontinuity is $s_{y}=\frac{2}{3}(\sqrt{a_{0}}-q_{0})$. Consequently,
leaving aside the growth of the extra leg on the top and bottom, the evolution
of the initial half solitons and the new stem becomes

$a(y,t)=\begin{cases}a_{0}&|y|>\frac{2}{3}(\sqrt{a_{0}}-q_{0})t\\\
(q_{0}+\sqrt{a_{0}})^{2}&|y|<\frac{2}{3}(\sqrt{a_{0}}-q_{0})t\end{cases},$
(51a)
$q(y,t)=\begin{cases}-\mathrm{sgn}(y)q_{0}&|y|>\frac{2}{3}(\sqrt{a_{0}}-q_{0})t\\\
0&|y|<\frac{2}{3}(\sqrt{a_{0}}-q_{0})t\end{cases}.$ (51b)

Comparing the modulation theory analytical solution (51) with parameter values
extracted from numerical simulation of the full KP equation (1) with
$\epsilon=1$, shown in figure 3, we see that (51) shows excellent agreement
with large $t$ numerical results in the bottom panels of figure 4.

Modulation theory can also analytically describe the growth of the new leg. We
focus on the upper left leg, with symmetrical results for the bottom left leg.
Since this leg must appear immediately to satisfy the jump conditions, we can
model its appearance as a partial soliton, a problem previously solved [61,
49, 50]. The partial soliton problem is a Riemann problem in the modulation
variables

$a_{\rm p}(y,0)=\begin{cases}a_{\rm n}&y<0\\\ 0&y>0\end{cases},\qquad q_{\rm
p}(y,0)=\begin{cases}q_{\rm n}&y<0\\\ q_{*}&y>0\end{cases},$ (52)

where $q_{*}$ can be determined by conservation of a Riemann invariant. The
solution to (52) is a simple wave with $R_{+}$ constant, and we can write down
an explicit solution for all values above the expanding stem where $a(y,t)$
and $q(y,t)$ are multivalued, i.e. for $y>\frac{2}{3}(\sqrt{a_{0}}-q_{0})t$,

$\displaystyle\sqrt{a_{\rm p}(y,t)}$
$\displaystyle=\begin{cases}0&\qquad\qquad\quad\,2q_{*}t<y\\\
\,\frac{3}{8}\left(2q_{*}-\frac{y}{t}\right)&\left(2q_{\rm
n}-\frac{2}{3}\sqrt{a_{\rm n}}\right)t<y<2q_{*}t\\\ a_{\rm
n}&\quad\frac{2}{3}(\sqrt{a_{0}}-q_{0})t<y<\left(2q_{\rm
n}-\frac{2}{3}\sqrt{a_{\rm n}}\right)t\end{cases},$ (53a) $\displaystyle
q_{\rm p}(y,t)$ $\displaystyle=q_{\rm n}+\sqrt{a_{\rm n}}-\sqrt{a_{\rm
p}(y,t)}.$ (53b)

Here $q_{*}=R_{+}=q_{\rm n}+\sqrt{a_{\rm n}}$ and represents the limit of
$q(y,t)$ for the $a=0$ edge of the simple wave. Note that the new leg is well-
defined, since the bottom (soliton) edge of the simple wave travels with
velocity $2q_{\rm n}-\frac{2}{3}\sqrt{a_{\rm
n}}=2\sqrt{a_{0}}-\frac{2}{3}q_{0}>\frac{2}{3}(\sqrt{a_{0}}-q_{0})$, so the
upper left leg always grows faster than the discontinuity with the stem. The
solution (53) is shown to accurately capture the large $t$ numerics in the top
panels of figure 4. To our knowledge, the evolution (53) has not been obtained
in previous work for V-shape initial conditions (49).

The above _Mach reflection_ behaviour occurs for $\sqrt{a_{0}}\geq q_{0}$. The
maximum amplitude of the intermediate (stem) region is obtained when
$q_{0}=\sqrt{a_{0}}$, which gives $a_{i}=(2\sqrt{a_{0}})^{2}=4a_{0}$,
precisely the maximum four-fold amplification first predicted by Miles [29].
When $q_{0}>\sqrt{a_{0}}$, the solution (51) breaks down, since the stem no
longer grows. We consider this case next.

###  5.2 Regular reflection

Figure 5: Numerical evolution of the KP equation (1) with $\epsilon=1$ for
V-shape initial conditions $\eqref{eq:ArrowheadIC}$ projected onto (50) on
$t\in\\{0,40,100\\}$ for $a_{0}=1$ and $q_{0}=1.4$. The resonance condition is
not met, so Mach reflection does not occur.

When $q_{0}>\sqrt{a_{0}}$ for (49), we have _regular reflection_. It is well
known that in this regime an X-shaped solution arises where two solitons
intersect with a small phase shift [15], which is higher-order [29] and
therefore not captured by modulation theory. A simulation of regular
reflection is shown in figure 5. The growth of the new left legs of the X can
be modelled identically as before in (53), with the only difference that
$a_{\rm n}=a_{0}$ and $q_{\rm n}=-q_{0}$, and the lower bound of the new
soliton is $y=0$. Consequently, $q_{*}=q_{0}+\sqrt{a_{0}}$. The partial
soliton solution (53) with the new parameters is favourably compared with
numerical simulation in figure 6.

The resulting X-shaped soliton, often called the O-type (meaning “ordinary”
[53]) soliton, is a travelling wave exact solution of the KP equation (1) with
velocity $(\overline{s_{x}},\overline{s_{y}})$. By the argument in section 4,
this structure as a whole should satisfy the jump conditions as given in
section 4. For $y>\overline{s_{y}}t$, we have two solitons with parameters
$(a_{0},\pm q_{0})$ while for $y<\overline{s_{y}}t$, we also have two solitons
with parameters $(a_{0},\mp q_{0})$. As a result, all $(n,m)=(2,2)$ modified
Rankine-Hugoniot jump conditions are trivially satisfied with no information
gained about $s_{y}$. However, we can still utilize jump conditions for $q$
from (15) to determine the velocity $\overline{s_{y}}$. Taking any two legs
with opposite signs for $q$, we obtain

$\overline{s_{y}}=\frac{\llbracket a/3+q^{2}\rrbracket}{\llbracket
q\rrbracket}=\frac{a_{0}/3+q_{0}^{2}-(a_{0}/3+q_{0}^{2})}{q_{0}-(-q_{0})}=0,$

as expected from symmetry.

Figure 6: Accuracy of the modulation prediction for the simulation from figure
5. The left two panels show contour plots for $t=30$ and $t=50$, with the
predictions for $\theta$ from (53) and (50) overlayed (dashed). The right
panel shows the accuracy of the amplitude prediction (dashed) compared to
numerical simulation (solid) for $t\in\\{40,80,120,160\\}$. The long-time
asymptotic prediction (dash-dotted) is also shown. The modulated prediction is
shown to be accurate, although the peak amplitude is actually $\approx 2.3$,
whereas superposition of modulation theory solutions predicts $a(0,t)=2$.

##  6 Discussion and conclusion

The main finding of this paper is that the KP resonant Y soliton can be
modelled as a shock solution to the soliton modulation equations. Remarkably,
the velocity of the Y soliton satisfies modified Rankine-Hugoniot jump
conditions for an infinite family of conservation laws. The justification for
our result is the fact that the modulation equations (22) represent the
soliton parameter evolution in the zero dispersion limit of the full KP
equation (1). In this limit, the rapid transition from the single, larger stem
of the resonant soliton to the two smaller legs limits to a single point of
discontinuity. This zero dispersion limit approach for defining admissible
discontinuous weak solutions is analogous to the vanishing viscosity approach
of admissible classical shocks in conservation law theory [65]. Our work is
the second known example of Whitham shocks [51], discontinuities in the
modulation variables that correspond to travelling wave solutions of the full
equation.

Since the Y soliton parameters are multivalued in the modulation equations,
and the fully two-dimensional modulation system reduces to an equivalent one-
dimensional system, we must introduce modified Rankine-Hugoniot conditions.
Namely, the fluxes and densities of the two smaller legs are added when
calculating the shock velocities. This is justified by the fact that the two
smaller legs summed together satisfy integral conservation laws when one
integrates over $x$.

Through Whitham modulation theory, we were able to analytically describe the
evolution of V-shape initial conditions (49). Unlike prior work, our approach
also determines the dynamic growth of the reflected legs as partial solitons.
These predictions were quantitatively verified to high accuracy in the large
$t$ regime by numerical simulation for the KP equation (1). In addition, we
identified the transition between Mach reflection and regular reflection for a
soliton incident upon an inward, compressive corner, with Mach reflection
occurring when $|q_{0}|<\sqrt{a_{0}}$.

Mach reflection is analogous to the phenomenon of Mach expansion, identified
in [49], where a soliton incident on an _outward_ corner creates a new,
straight _lower-amplitude_ stem whenever $|q_{0}|<\sqrt{a_{0}}$. Within the
Mach reflection and expansion regimes, a new line soliton is formed
perpendicular to the wall. Outside of these regimes, the interaction of a
soliton with an outward or inward corner is very different, with the initial
partial solitons evolving independently from each other.

The results in this paper are (to our knowledge) the first application of
modulation theory to a fully two-dimensional evolution, thereby opening new
doors for future work studying soliton resonance in other two-dimensional
equations. Modulation theory does not depend on integrability or, relatedly,
the existence of a large family of multi-soliton solutions. Consequently, this
work can be generalized to any system with a known line soliton solution. For
example, the two-dimensional Benjamin-Ono (2DBO) equation, which models
internal water waves [66], has numerically-observed resonance-like behaviour
[37] and line soliton solutions for which modulation equations have been
calculated [67]. Few analytical tools and no known exact multi-soliton
solutions exist for the 2DBO equation, but the generalization of our approach
seems like a promising avenue for research in this area.

## Acknowledgments

Authors thank Jonathan Hauenstein for helpful discussions concerning
polynomial systems. The work of MAH and SR was supported by NSF grant
DMS-1816934. The work of GB was supported by NSF grant DMS-2009487.

##  Appendix A Reducibility of multidimensional modulations

In this appendix we prove Theorem 3.1. Namely, we will show that when $a$ and
$q$ describe a modulated KP soliton (2), the system (12) can be reduced to a
(1+1)-dimensional system, assuming differentiable modulations.

###### Proof.

The full modulation system (12) is equivalent to (22) under the assumption
that $a_{x}\equiv q_{x}\equiv 0$, which followed from the ansatz (14). This
assumption can be relaxed by using a more general line soliton ansatz as

$u_{0}(\theta,x,y,t;\epsilon)=a(x,y,t)\operatorname{sech}^{2}\left(\frac{1}{\psi}\sqrt{\frac{a(x,y,t)}{12}}\frac{\theta}{\epsilon}\right),$
(54)

where $\psi=\psi(a,q)$ is nonzero, real-valued, and differentiable. If the
soliton is unmodulated (i.e. $a(x,y,t)\equiv a$ and $q(x,y,t)\equiv q$), (54)
must be consistent with (2). Since the argument of $\operatorname{sech}^{2}$
in (54) has been divided by $\psi$, it must also be multiplied by $\psi$. We
incorporate $\psi$ into the derivatives of the fast phase $\theta/\epsilon$ so
that

$\theta_{x}=\psi(a,q),\qquad\theta_{y}=q\psi(a,q),\qquad\theta_{t}=-\psi(a,q)\left(\frac{a}{3}+q^{2}\right),$
(55)

where $a=a(x,y,t)$ and $q=q(x,y,t)$ are the modulation parameters. For the
ansatz in this work (14), we selected $\psi=1$, while for a previous work
focusing on $y$-independent modulations we selected $\psi=1/q$ [50]. The
compatibility conditions $\theta_{xy}=\theta_{yx}$ and
$\theta_{xt}=\theta_{tx}$ for (55) yield, respectively

$\displaystyle\psi_{a}a_{y}+\psi_{q}q_{y}$
$\displaystyle=q\psi_{a}a_{x}+(q\psi_{q}+\psi)q_{x},$ (56a)
$\displaystyle\psi_{a}a_{t}+\psi_{q}q_{t}$
$\displaystyle=-\left(c\psi_{a}+\frac{\psi}{3}\right)a_{x}-\left(c\psi_{q}+2q\psi\right)q_{x}.$
(56b)

We require consistency with the full modulation system (12), so we insert the
equations for $a_{t}$ and $q_{t}$ from (12) into (56b) to obtain an equation
in which $(a,q)$ only depend on the spatial variables $(x,y)$. Using
simplifying calculations together with (56a), we obtain the system

$\begin{bmatrix}\psi_{a}&\psi_{q}\\\
\psi_{q}&4a\psi_{a}\end{bmatrix}\begin{bmatrix}a_{y}\\\
q_{y}\end{bmatrix}=q\begin{bmatrix}\psi_{a}&\psi_{q}\\\
\psi_{q}&4a\psi_{a}\end{bmatrix}\begin{bmatrix}a_{x}\\\
q_{x}\end{bmatrix}+\begin{bmatrix}0&\psi\\\
\psi&0\end{bmatrix}\begin{bmatrix}a_{x}\\\ q_{x}\end{bmatrix}\,.$ (57)

We can always solve for either $(a_{x},q_{x})$ or $(a_{y},q_{y})$ in equation
(57), since for any $\psi$ either the left hand side or the right hand side
matrix in (57) is invertible. We show this by assuming both determinants are
zero and generating a contradiction. Let us rewrite (57) as
$A\bm{v}_{y}=B\bm{v}_{x}$, where $\bm{v}=[a\;q]^{T}$ and $A,B$ are the
appropriate matrices from (57). We assume that

$\displaystyle\det{A}$ $\displaystyle=4a\psi_{a}^{2}-\psi_{q}^{2}=0,$ (58a)
$\displaystyle\det{B}$
$\displaystyle=4a\psi_{a}^{2}-\psi_{q}^{2}-\psi^{2}-2q\psi_{q}\psi=-\psi^{2}-2q\psi_{q}\psi=0.$
(58b)

The condition (58b) for nonzero $\psi$ implies that $\psi+2q\psi_{q}=0$. This
ordinary differential equation is solved as $\psi=R(a)/\sqrt{q}$ for any
$R(a)$. But inserting this into (58a), we obtain an equation for $R$ that is
satisfied when $4q\sqrt{a}R^{\prime}(a)=\pm R(a)$. This contradicts the fact
that $R=R(a)$ is only a function of $a$, so we cannot find a solution for
$\psi$ that satisfies (58) and sets both determinants to zero. Therefore, for
any $\psi$ either $A$ or $B$ must have a nonzero determinant and therefore be
invertible. ∎

Consequently, in (57) we can solve for either $(a_{x},q_{x})$ or
$(a_{y},q_{y})$. We can then locally substitute whichever variables were
solved for into the full modulation system (12), reducing (12) to a
(1+1)-dimensional system.

## References

  * [1] J. Scott Russell.  Report on waves.  In Report of the fourteenth meeting of the British Association for the Advancement of Science, pages 311–390. John Murray, London, September 1845.
  * [2] E. Mach and J. Wosyka.  Über einige mechanische Wirkungen des electrischen Funkens.  Sitzungsber. Akad. Wiss. Wien, 72:44–52, 1875.
  * [3] P. Krehl and M. van der Geest.  The discovery of the Mach reflexion effect and its demonstration in an auditorium.  Shock Waves, 1:3–15, 1991.
  * [4] B. B. Kadomtsev and V. I. Petviashvili.  On the stability of solitary waves in weakly dispersing media.  Sov. Phys. Dokl., 15(6):753–756, 1970.
  * [5] M. J. Ablowitz and P. A. Clarkson.  Solitons, Nonlinear Evolution Equations and Inverse Scattering.  Cambridge University Press, 1991.
  * [6] M. Ablowitz and D. Baldwin.  Nonlinear shallow ocean wave soliton interactions on flat beaches.  Physical Review E, 86:036305, 2012.
  * [7] M.J. Ablowitz and H. Segur.  On the evolution of packets of water waves.  J. Fluid Mech., 92(04):691–715, 1979.
  * [8] R. Grimshaw.  Evolution equations for long, nonlinear internal waves in stratified shear flows.  Stud. Appl. Math., 65(2):159–188, 1981.
  * [9] C. Yuan, R. Grimshaw, E. Johnson, and X. Chen.  The propagation of internal solitary waves over variable topography in a horizontally two-dimensional framework.  Journal of Physical Oceanography, 48:283–300, February 2018.
  * [10] E. Infeld and G. Rowlands.  Nonlinear Waves, Solitons, and Chaos.  Cambridge UP, 2001.
  * [11] M. Funakoshi.  Reflection of obliquely incident solitary waves.  J. Phys. Soc. Japan, 49(6):2371–2379, 1980.
  * [12] Mitsuhiro Tanaka.  Mach reflection of a large-amplitude solitary wave.  J. Fluid Mech., 248:637–661, 1993.
  * [13] Hidekazu Tsuji and Masayuki Oikawa.  Oblique interaction of solitons in an extended Kadomtsev-Petviashvili equation.  J. Phys. Soc. Japan, 76(8):084401, 2007.
  * [14] G. Biondini, K.-I. Maruno, M. Oikawa, and H. Tsuji.  Soliton interactions of the Kadomtsev-Petviashvili equation and generation of large-amplitude water waves.  Stud. in Appl. Math., 122(4):377–394, 2009.
  * [15] Y Kodama, M Oikawa, and H Tsuji.  Soliton solutions of the KP equation with V-shape initial waves.  Jour. of Phys. A: Math. and Theo., 42(31):312001, July 2009.
  * [16] Wenwen Li, Harry Yeh, and Yuji Kodama.  On the Mach reflection of a solitary wave: revisited.  J. Fluid Mech., 672:326–357, 2011.
  * [17] C.-Y. Kao and Y. Kodama.  Numerical study of the KP equation for non-periodic waves.  Mathematics and Computers in Simulation, 82(7):1185–1218, March 2012.
  * [18] S. Chakravarty, T. McDowell, and M. Osborne.  Numerical studies of the KP line-solitons.  Communications in Nonlinear Science and Numerical Simulation, 44:37–51, March 2017.
  * [19] Jeffrey Knowles and Harry Yeh.  Fourfold amplification of solitary-wave mach reflection at a vertical wall.  J. Fluid Mech., 861:517–523, February 2019.
  * [20] F. R. Gilmore, M. S. Plesset, and H. E. Crossley.  The analogy between hydraulic jumps in liquids and shock waves in gases.  J. of Appl. Phys., 21(3):243–249, 1950.
  * [21] Paul Henri Perroud.  The solitary wave reflection along a straight vertical wall at oblique incidence.  Institute of Engineering Research. Wave Research Laboratory. Technical report ser. 99, University of California, Berkeley, 1957\.
  * [22] L.F. Henderson and A. Lozzi.  Experiments on transition of Mach reflexion.  J. Fluid Mech., 68:139–155, 1975.
  * [23] W. K. Melville.  On the Mach reflexion of a solitary wave.  J. Fluid Mech., 98(2):285–297, 1980.
  * [24] M. Tanaka.  Mach reflection of a large-amplitude solitary wave.  J. Fluid Mech., 248:637–661, 1993.
  * [25] H. Yeh, W. Li, and Y. Kodama.  Mach reflection and kp solitons in shallow water.  Eur. Phys. J. Spec. Top., 185:97–111, 2010.
  * [26] Y. Kodama and H. Yeh.  The KP theory and Mach reflection.  J. Fluid Mech., 800:766–786, 2016.
  * [27] C. Wang and R. Pawlowicz.  Oblique wave-wave interactions of nonlinear near-surface internal waves in the Strait of Georgia.  J. of Geophys. Res.: Oceans, 117(C6), 2012.
  * [28] K. Shimizu and K. Nakamaya.  Effects of topography and Earth’s rotation on the oblique interaction of internal solitary-like waves in the Andaman Sea.  Journal of Geophysical Research: Oceans, 122(9):7449–7465, 2017\.
  * [29] J. Miles.  Obliquely interacting solitary waves.  J. Fluid Mech., 79:157–169, 1977.
  * [30] J. Miles.  Resonantly interacting solitary waves.  J. Fluid Mech., 79:171–179, 1977.
  * [31] Ryogo Hirota.  The Direct Method in Soliton Theory.  Cambridge University Press, Cambridge, 2004.
  * [32] Y. Kodama.  Young diagrams and $N$-soliton solutions of the KP equation.  J. Phys. A, 37:11169–11190, 2004.
  * [33] G. Biondini.  Line soliton interactions of the Kadomtsev-Petviashvili equation.  Phys. Rev. Lett., 99:064103, August 2007.
  * [34] S. Chakravarty and Y. Kodama.  Soliton solutions of the KP equation and application to shallow water waves.  Stud. Appl. Math., 123(1):83–151, 2009.
  * [35] Yuji Kodama.  KP solitons in shallow water.  J. Phys. A, 43(43):434004, October 2010.  arXiv: 1004.4607.
  * [36] Jen-Hsu Chang.  Soliton interaction in the modified Kadomtsev-Petviashvili-(II) equation.  Applicable Analysis, 98:2589–2603, 2019.
  * [37] H. Tsuji and M. Oikawa.  Oblique interaction of internal solitary waves in a two-layer fluid of infinite depth.  Fluid Dynamics Research, 29(29):251–267, 2001.
  * [38] M. Oikawa and H. Tsuji.  Oblique interactions of weakly nonlinear long waves in dispersive systems.  Fluid Dynamics Research, 38(12):868–898, 2006.
  * [39] Y. Matsuno.  Oblique interaction of interfacial solitary waves in a two-layer deep fluid.  Proc. Roy. Soc. A, 454(1971):835–856, 1998.
  * [40] G.B. Whitham.  Linear and Nonlinear Waves.  Wiley Interscience, 1974.
  * [41] G.A. El and M.A. Hoefer.  Dispersive shock waves and modulation theory.  Physica D, 333:11–65, October 2016.
  * [42] Peter D. Lax and C. David Levermore.  The small dispersion limit of the korteweg-de vries equation. i.  Communications on Pure and Applied Mathematics, 36(3):253–290, 1983\.
  * [43] Peter D. Lax and C. David Levermore.  The small dispersion limit of the korteweg-de vries equation. ii.  Communications on Pure and Applied Mathematics, 36(5):571–593, 1983\.
  * [44] Peter D. Lax and C. David Levermore.  The small dispersion limit of the korteweg-de vries equation. iii.  Communications on Pure and Applied Mathematics, 36(6):809–829, 1983\.
  * [45] M.J. Ablowitz, G. Biondini, and Q. Wang.  Whitham modulation theory for the Kadomtsev–Petviashvili equation.  Proc. R. Soc. A, 473(2204):20160695, August 2017.
  * [46] Mark J. Ablowitz, Gino Biondini, and Igor Rumanov.  Whitham modulation theory for (2+1)-dimensional equations of Kadomtsev-Petviashvili type.  J. Phys. A, 51(21):215501, April 2018.
  * [47] T. Grava, C. Klein, and G. Pitton.  Numerical study of the Kadomtsev–Petviashvili equation and dispersive shock waves.  Proc. R. Soc. A, 474(2210):20170458, February 2018.
  * [48] G. Biondini, M.A. Hoefer, and A. Moro.  Integrability, exact reductions and special solutions of the KP–Whitham equations.  Nonlinearity, 33:4114–4132, 08 2020.
  * [49] Samuel Ryskamp, Michelle D. Maiden, Gino Biondini, and Mark A. Hoefer.  Evolution of truncated and bent gravity wave solitons: the Mach expansion problem.  J. Fluid Mech., 909:A24, 2021.
  * [50] Samuel J. Ryskamp, Mark A. Hoefer, and Gino Biondini.  Oblique interactions between solitons and mean flows in the Kadomtsev-Petviashvili equation, 2020.
  * [51] P. Sprenger and Mark A. Hoefer.  Discontinuous shock solutions of the whitham modulation equations as zero dispersion limits of traveling waves.  Nonlinearity, 33:3268, 2020.
  * [52] Sarbarish Chakravarty and Yuji Kodama.  Construction of KP solitons from wave patterns.  Journal of Physics A: Mathematical and Theoretical, 47(2):025201, dec 2013.
  * [53] Y. Kodama.  Solitons in Two-Dimensional Shallow Water.  SIAM, Philadelphia, PA, 2018.
  * [54] G. Biondini and S. Chakravarty.  Soliton solutions of the Kadomtsev-Petviashvili ii equation.  J. Math. Phys., 47(3):033514, March 2006.
  * [55] G. Biondini and Y. Kodama.  On a family of solutions of the Kadomtsev–Petviashvili equation which also satisfy the toda lattice hierarchy.  J. Phys. A, 36(42):10519–10536, October 2003.
  * [56] Harvey Segur and Allan Finkel.  An analytical model of periodic waves in shallow water.  Studies in Applied Mathematics, 73(3):183–220, 1985.
  * [57] V.A. Mironov, A.I. Smirnov, and L.A. Smirnov.  Dynamics of vortex structure formation during the evolution of modulation instability of dark solitons.  J. Exp. Theor. Phys., 112:46–59, 2011.
  * [58] G.B. Whitham.  Non-linear dispersive waves.  Proc. Roy. Soc. Ser. A, 283(1393):238–261, 1965.
  * [59] S.C. Anco, M.L. Ganadrias, and E. Recio.  Conservation laws, symmetries, and line soliton solutions of generalized KP and Boussinesq equations with $p$-power nonlinearities in two dimensions.  Theoretical and Mathematical Physics, 197(1):1393–1411, 2018.
  * [60] Stephen C. Anco and Elena Recio.  Topological charges and conservation laws involving an arbitrary function of time for dynamical PDEs.  Proc. Roy. Soc. A, 477:2020442, 2021.
  * [61] J.C. Neu.  Singular Perturbation in the Physical Sciences.  Am. Math. Soc., Providence, 2015.
  * [62] D.E. Baldwin and W. Hereman.  A symbolic algorithm for computing recursion operators of nonlinear partial differential equations.  Int’l J. of Comp. Math., 87(5):1094–1119, 2010.
  * [63] L. Ostrovsky and Y.A. Stepanyants.  Kinematics of interacting solitons in two-dimensional space.  Russ. J. of Earth Sci., 20(ES4007), Jul 2020.
  * [64] Adam Mahdi, Claudio Pessoa, and Jonathan D. Hauenstein.  A hybrid symbolic-numerical approach to the center-focus problem.  Journal of Symbolic Computation, 82:57–73, 2017.
  * [65] P. LeFloch.  Hyperbolic systems of conservation laws.  Springer, Boston, 2002.
  * [66] M. J. Ablowitz and H. Segur.  Long internal waves in fluids of great depth.  Stud. Appl. Math., 62(3):249–262, 1980.
  * [67] Mark Ablowitz, Gino Biondini, and Qiao Wang.  Whitham modulation theory for the two-dimensional Benjamin-Ono equation.  Physical Review E, 96(3), Sep 2017.