Librating Kozai-Lidov Cycles with a Precessing Quadrupole Potential are Analytically Approximately Solved

The very long-term evolution of the hierarchical restricted three-body problem with a slightly aligned precessing quadrupole potential is investigated analytically for librating Kozai-Lidov cycles (KLCs). \citet{klein2023} presented an analytic solution for the approximate dynamics on a very long timescale developed in the neighborhood of the KLCs fixed point where the eccentricity vector is close to unity and aligned (or anti aligned) with the quadrupole axis and for a precession rate equal to the angular frequency of the secular Kozai-Lidov Equations around this fixed point. In this Letter, we generalize the analytic solution to encompass a wider range of precession rates. We show that the analytic solution approximately describes the quantitative dynamics for systems with librating KLCs for a wide range of initial conditions, including values that are far from the fixed point which is somewhat unexpected. In particular, using the analytic solution we map the strikingly rich structures that arise for precession rates similar to the Kozai-Lidov timescale (ratio of a few).


INTRODUCTION
In this Letter we generalize and study the applicability of the analytical solution introduced by Klein & Katz (2023) for the long-term evolution of a test particle in a Keplerian orbit perturbed by a slightly aligned precessing quadrupole potential.The system under consideration involves a test particle orbiting a central mass M on a Keplerian orbit with semimajor axis a that is perturbed by an external quadrupole potential given by where Φ 0 is constant and ĵouter is a time dependent unit vector which precesses around the z -axis at a constant rate β with a constant inclination α: where τ ≡ t tsec and t sec = For β = 0 the potential remains constant over time and the system represents the scenario of a hierarchical three-body system, where the gravitational potential has been expanded to quadrupole order in terms of the separations ratio and averaged over the outer period.This problem has been analytically solved using the secular approximation (Kozai (1962); Lidov (1962), for a recent review see Naoz (2016)) resulting with the periodic Kozai-Lidov cycles (KLCs) which reach high eccentricities for high initial inclinations 1 .
The emergence of high eccentricity for the inner binary of a hierarchical system in a KLC has been proposed as a crucial factor in the possible formation channel of a wide variety of astrophysical phenomena such as hot/warm planets (Naoz et al. 2011;Katz et al. 2011;Fabrycky & Tremaine 2007), planets around white dwarfs (Muñoz & Petrovich 2020;O'Connor et al. 2020;Stephan et al. 2021), Type Ia supernovae through the merger or collision of white dwarfs in multiple systems (Thompson 2011;Katz & Dong 2012), gravitational wave emission through the merger of black holes or neutron stars (An-tonini & Perets 2012;Liu & Lai 2018a) and the formation of close binaries (Antonini & Perets 2012).
Previous numerical investigations have demonstrated that in systems with higher multiplicities, whether they include an additional perturber, account for the effect of the host galaxy, or consider the asphericity of star clusters, the phase space that allows for high eccentricities is more extensive when compared to equivalent triple systems (Pejcha et al. 2013;Hamers et al. 2015;Hamers 2017;Petrovich & Antonini 2017;Fang et al. 2018;Grishin et al. 2017;Liu & Lai 2018b;Safarzadeh et al. 2019;Hamers & Safarzadeh 2020;Bub & Petrovich 2020;Grishin & Perets 2022).The simplest case for rich dynamics in a higher multiplicity system is a case where the inner binary's angular momentum is negligible (the test particle assumption) and there are no KLC dynamics for the other binaries in the systems.This simplest case of higher multiplicity (with some restricting assumptions) is equivalent to the precessing quadrupole studied in this Letter, Equations 1 and 2, and numerical explorations have validated that even for low alignments it exhibits a wide range of initial conditions leading to high eccentricity outcomes (Hamers & Lai 2017;Petrovich & Antonini 2017).
The dynamics of the test particle can be parameterized by two dimensionless orthogonal vectors j = J/ √ GM a, where J is the specific angular momentum vector, and e a vector pointing in the direction of the pericenter with magnitude e.In the secular approximation, where the equations are averaged over the orbit, a is constant with time while j and e evolve according to As mentioned above, in the periodic analytically solved KLCs the external potential is constant in time implying an axisymmetric potential and admitting a constant of the motion, j • ĵouter .In this case a second constant of motion can be obtained from the double averaged potential and j • ĵouter The dynamics of a KLC are determined by the values of the constants C K and j • ĵouter .When C K < 0, the argument of pericenter of the Keplerian orbit librates around π 2 or − π 2 (librating cycles), and when C K > 0, it goes through all values [0, 2π] (rotating cycles).
In the problem we study here with β > 0, i.e ĵouter is precessing at a constant rate, C K and j • ĵouter are no longer constant.We note that in the presence of the third order term (octupole) in the expansion of the potential a change in C K and j• ĵouter leading to high eccentricities and orbit flips can be obtained (Naoz et al. 2011;Katz et al. 2011;Lithwick & Naoz 2011;Teyssandier et al. 2013;Li et al. 2014;Naoz 2016;Tan et al. 2020;Lei 2022;Huang & Ji 2022).In the precessing potential studied here the change in C K and j • ĵouter is obtained at the quadrupole level (Petrovich & Antonini 2017).
A global constant of motion of Equations 2-3 can be derived as follows (equivalent to the constant Hrot presented in Petrovich & Antonini (2017)).Time independent equations can be obtained by transforming to a reference frame rotating with the perturbing potential.This is equivalent to adding (−βj z ) to the Hamiltonian (Tremaine & Yavetz 2014) and implies the following constant of motion (equal to the resulting time independent Hamiltonian up to additive and multiplicative constants) Since this expression is independent of the reference frame it represents a global constant of Equations 2-3 as can be readily checked directly2 .This constant applies for all initial conditions and any value of α.
As outlined in Klein & Katz (2023), when 0 < α ≪ 1 the dynamics on short time scales adheres to the conventional KLC with two conserved quantities: j • ĵouter ≈ j z and C K (Equation 4).On longer time scales j z and C K evolve (satisfying the global constant of Equation 5).In this Letter, our emphasis is directed towards the circumstances conducive to high eccentricities, which correspond to low |j z | (extremely high eccentricities correspond to a zero crossing of j z , e.g Naoz et al. (2011)).A measure of the range of initial j z that can lead to a zero crossing outcome is the magnitude of variation observed in j z in the vicinity of j z = 0 (Katz et al. 2011;Luo et al. 2016;Haim & Katz 2018).This metric will serve as the primary quantitative parameter that we will examine, both numerically and analytically.
In Klein & Katz (2023) we presented an analytic solution to the KLC averaged equations of this problem for α ≪ 1 in a restricted regime close to the KLC fixed point of 2 and j = 0 (i.e e = 1).The effective averaged equations and analytic solution in Klein & Katz (2023) were derived for a specific value of precession rate which is the natural resonant frequency of the equations of motion in the neighborhood of this fixed point (e close to unity and aligned or anti-aligned with ĵouter and j ≈ 0).The success of the solution to capture and approximately describe the long-term evolution of the dynamics was shown for examples with β close to β 0 and initial conditions close to the fixed point.
In this Letter, we generalize the KLC averaged equations and their analytic solution for a wider range of β values and compare it to the numerical results for a wide range of α and β.We explore a large phase space of initial conditions close and far from the fixed point C K = −1.5 but we restrict our analysis to librating KLCs (i.e C K < 0 at all times), α ≪ 1 and |j z | ≪ 1.

NUMERICAL RESULTS
We numerically integrate Equations 2-3 (up to τ = 300), with varying initial conditions for the vectors e and j, as well as different values of α and β.Our analysis is restricted to simulations where the parameter C K remains negative throughout its temporal evolution, implying persistent libration in the KLC conditions.With an emphasis on the low |j z | region (favorable for inducing high eccentricity), we set the initial value of j z to zero.This implies that in all the results shown, the initial value of C K , C 0 K , is equal to the global constant in Equation 5 up to first order in α.To cover the multidimensional phase space of initial conditions, for each combination of α and β values under investigation, we generate 1000 instances of initial conditions.Every instance is produced by randomly selecting three values for e x , e y , and j x from a uniform distribution within the interval [−1, 1].We then determine e z and j y by demanding the conditions e • j = 0 and e 2 + j 2 = 1, selecting one solution at random from the two possible solutions.If no solution exists -a new random selection of e x , e y , and j x is made instead, implying 1000 instances of viable initial conditions.If C 0 K > 0 or C K becomes positive during the evolution we discard the simulation from the analysis.The number of presented integrations is therefore less than 1000 and is indicated explicitly for each result shown (top of the relevant subplots in the figures below).
The amplitude of change in j z throughout the numerical integration (∆j z = j max Note the following striking features of the numerical results shown in Figure 1: (a) j z can change considerably (even) for low α, an effect that is quickly suppressed when β becomes larger than β 0 .(b) ∆j z is determined mainly by C 0 K .(c) The relationship between ∆j z and C 0 K shows intricate and non-uniform patterns: Various correlations are observed, distinct subdivisions into separate sub-categories are evident, notable levels of dispersion occur, there are instances of remarkably narrow correlation, and abrupt localized increments emerge.
The green circles in Figure 1 denote results from the analytic solution outlined in the upcoming sections.Note the strong resemblance for most (but not all) of the distinctive characteristics between the two solutions.A detailed comparison is discussed in section 4.

AVERAGED EQUATIONS AND ANALYTIC SOLUTION
As demonstrated in Klein & Katz (2023), close to the fixed point C K = −1.5, up to first order in α and within the low |j z | regime, Equations 2-3 turn out to have a structure of a sinusoidal-driven harmonic oscillator equation for e x and e y , under the assumption of a constant e z .This simplification leads to the following equations: ëx = ω 2 0 (L cos (ωτ ) − e x ) and ëy = ω 2 0 (L sin (ωτ ) − e y ) where ω = β, L = e z α and ω 0 = β 0 |e z | (see Equations 4-9 in Klein & Katz (2023)).
As a consequence, a measure of the proximity to resonance for these driven harmonic oscillators involves evaluating the difference ω 2 0 − ω 2 .In the context of a standard driven harmonic oscillator, both ω 0 and ω remain constant over time, thereby establishing a fixed state of proximity to resonant behavior.However, in our examined scenario, ω = β remains unchanged while ω 0 = β 0 |e z | can slowly evolve, in addition to variations in the instantaneous e z induced by the KLC itself.In the following we shall average the equations over KLCs.The deviation from resonance, ω 2 0 −ω 2 , is represented by one of the slow dynamical variables that we follow (see Equation 9 below).Non-trivial dynamics are obtained when the system approaches or even crosses resonance as the KLC-averaged value of e z gradually changes over time.Klein & Katz (2023))

Ansatz (as in
For convenience we lay out in this subsection definitions and equations that were given in Klein & Katz (2023): We use the following ansatz for the vector e in where ϕ is a slowly dynamically evolving phase and We use the following dynamical parameter measuring the proximity to resonance in the driven harmonic oscillator (with ω = β fixed but ω 0 = β 0 |ē z | slowly evolving with time), where ēz is the averaged value of e z over KLC which satisfies

Averaged Equations
In what follows the equations presented are new, with β not restricted to equal β 0 but assumed to be in the vicinity of β 0 |e z |.
Using the following slow variables where the following set of ODEs is obtained: where ˙denotes a derivative with respect to τ .Note that a + b is a constant of motion and the equations for δ, s and c form a closed subset.
Using the demand that j • e = 0 we deduce that the following combination of j z and δ is constant: enabling j z to be readily obtained at any time utilizing the initial conditions and the value of δ.

Analytic Solution
The averaged equations, Equations 16-18, admit two constants of the motion, denoted C 1 and C 2 , implying that the evolution of the three variables δ, a−b and θ − ϕ is periodic.
The resulting evolution of δ is equivalent to the dynamics of a particle moving in a one dimensional potential with a constant energy where The pair of constants, C 1 and C 2 , which determine V and E offer a comprehensive analytical solution for the evolution of δ.Consequently, this solution allows for the derivation of the evolution of all the slowly evolving variables.Note that the energy and potential are defined slightly different than in Klein & Katz (2023) (Equations 29-30 there) with an extra added constant 1 2 β0 β C 1 2 so E would depend on C 2 only.This addition of a constant value does not alter the results which are presented as a function of C 1 and C 2 and could have been added there.
Equations 22 and 23 determine the amplitude of change in δ analytically.Using Equation 19 we obtain the analytic prediction for ∆j z which is plotted using green circles in Figure 1.
We note that at first glance, the analytic solution remarkably agrees with the general trends and most of the striking features and patterns identified in the numerical outcomes outlined in section 2. It effectively captures the principal observations solely using the initial conditions: (a) The large j z changes are approximately quantitatively reproduced for a large phase space of C 0 K .In particular, the suppression of the effect as β becomes larger than β 0 is recovered.(b) Like in the numerical results ∆j z is determined mainly by C 0 K .(c) The rich variety of types of dependence of ∆j z on C 0 K is mainly generated also by the analytic solution: The distinct division to separate sub-classes and the regions where a significant scatter is observed.Note that there are features that are not reproduced.In particular, a strikingly narrow correlation is observed in the numerical results for C 0 K > −0.6 in some of the panels and is not reproduced.This very tight correlation in the numerical results is surprising and warrants further investigation which is beyond the scope of this work.Additionally, in the two lower left panels there are localized increases (at C 0 K ∼ −0.3) that are not recovered.

Distinct Qualitative Pairs of V and E
The analytic solution assigns each set of parameters and initial conditions to one of four distinct classes defined below.The potential V in Equation 23has two distinct shapes depending on whether C 1 is smaller or larger than a critical value (using Equation 13) the potential V has no maxima and one minimum.If C 1 > C crit 1 the potential has a local maximum (and two minima) and we denote by V max the value of the potential at this local maximum and by δ Vmax the δ at which V max is obtained.The properties of Equation 23imply that δ Vmax > 0.
The phase space of initial conditions is separated to four distinct classes (colors correspond to Figures 2, 3  and 5): , E < V max , δ > δ Vmax : see example in the right lower panel of Figure 2.
An immediate relation between these classes of initial conditions and ∆j z seems to emerge from the examples shown in Figure 2: Large changes in j z which correspond to large changes in δ are expected for classes 1 and 2, and not for classes 3 and 4. As shown in section 4 this is true for a wide range of initial conditions and parameters.

COMPARISON OF TRENDS BETWEEN ANALYTIC AND NUMERICAL SOLUTIONS
In the following subsections we use the analytic solution to explain the main different structures shown in Figure 1.
The numerical results shown in Figure 1 are reproduced in Figure 3 with different colors corresponding to the analytical classification of the initial conditions defined in section 3.4.We note several striking observations: (a) Different numerical features correspond to distinct analytical classes (note a slight exception for Class 2 (blue) in some of the panels which will be discussed in section 4.3).(b) As claimed in section 3.4, large changes in j z mostly correspond to Class 1 (red) or Class 2 (blue).(c) The subset characterized by low values of ∆j z at the negative end of C 0 K predominantly aligns with Class 4 (cyan).(d) The region demonstrating a notably tight correlation appearing as a rightward tail (high C 0 K and low ∆j z ) corresponds to Class 3 (magenta).(e) Initial conditions identified as Class 1 (red) have typically the largest scatter in ∆j z .

α dependence
The right column of panels in Figures 1 and 3 together with the middle left panel show the results with a constant value of β = 2.4 and a logarithmically increasing α from 0.125 • to 4 • .
The effect of varying α is three fold: • As α increases -the potential amplitude of ∆j z increases.As can be seen in these subplots in Figure 1 the analytic solution captures quantitatively the dependence of the amplitude of the effect on α.Notably, the y-axis scale varies significantly across these subplots, and the green circles track the trends observed in the black crosses.
• Increasing α induces changes in the color distributions shown in the maps, with the portion of initial conditions attributed to Class 1 (red) increasing with α and reducing the number of occurrences of classes 3 (magenta) and 4 (cyan).The heightened potential amplitude of ∆j z with increasing α can be ascribed to a dual effect: the prevalence of Class 1 cases with a greater amplitude of δ, and the influence of the α 2 3 term in Equation 19, which translates changes in δ into ∆j z .
• As α increases the count N of librating cases experiences a decline.This phenomenon is rooted in the intensified perturbing potential, causing more vigorous fluctuations in C K and subsequently leading to an increased occurrence of cases where

β dependence
The left column of panels in Figures 1 and 3 show the results with a constant value of α = 1 • and a range of β values from 1.2 to 3.5 (including β 0 ≈ 2.9) for which considerable change in ∆j z (≳ 0. The change of β predominantly impacts the distribution of the different classes based on the initial conditions (i.e C 0 K ).As evidenced by these subplots in Figure 3, at β = 1.2 the majority of initial conditions fall into Class 4 (cyan).As β increases, the proportion of Class 4 (cyan) gradually diminishes, giving way to an increasing fraction of Class 1 (red) and 3 (magenta).Specifically, for β values of 1.8 and 2.4, a substantial portion of the phase space is occupied by Class 1 (red), enabling significant changes in δ and consequently yielding high values of ∆j z .When β ≥ β 0 no instances of Class 2 (blue) or 4 (cyan) are present, as actual initial conditions for e invariably satisfy 0 ≤ ē2 z ≤ 1.Consequently, from Equation 9, it becomes evident that δ < 0 for β ≥ β 0 , restricting the cases to only classes 1 (red) or 3 (magenta).Further elevating β away from β 0 , for instance, to β = 3.5, leads to all cases being categorized as Class 3 (magenta).This phenomenon emerges as a common trait of sufficiently large β, as an increased minimal value for δ 2 is observed with rising β (see Equation 9), and consequently, C 1 (linked to δ 2 , see Equation 20) quickly surpasses C crit 1 for all initial conditions, making Class 1 (defined by C 1 < C crit 1 ) unattainable.Generally, these subplots in Figure 1 demonstrate that the analytic solution captures quantitatively the dependence of ∆j z on β as the green circles follow the generic trends of the black crosses (although failing at low values of C 0 K for β = 1.2).The highest values attained by ∆j z in these panels is shown as a function of β in Figure 4 along with results from additional values of β for α = 1 • and the corresponding analytic results.As in Figure 1 black crosses denote numerical results and green circles denote the analytic prediction using the initial conditions.Note the significant elevation of maximal ∆j z in the range 1 ≲ β ≲ 3 and the steep decline for β > β 0 as mentioned above.The analytic solution approximately reproduces the numerical result.A large deviation is observed at β ∼ 1 and corresponds to the large values of ∆j z in the analytic prediction around C 0 K ∼ −0.3 as can be seen in the upper left panel of Figure 1 showing the results for β = 1.2.Note that for lower C 0 K values the agreement between the solutions is much better for this β.

Detailed comparison for a representative case -
A detailed comparison between the numerical and analytical solutions for a representative case of α = 1 • and β = 2.4 is given in Figure 5.We present the numerical results (left column) and analytical results (right column) with the aforementioned variable of ∆j z (top panels), and we add the values composing it: j max z (middle panels) and j min z (bottom panels).This set of α and β values is chosen since high ∆j z is reached and rich structures are obtained as can be seen in the middle left panel of Figure 3.
As observed in the top right panel, Class 4 cases are a distinct subcategory with a low ∆j z .As apparent in the middle and bottom panels -this distinction appears .This can be explained by the analytic solution as follows.First, note that using Equation 19, j max z (j min z ) corresponds to δ min (δ max ).Second, note that when transitioning from Class 4 to Classes 1 or 2, there is a discontinuity in δ min but not in δ max due to the shape of the potential as can be seen in Figure 2.
A prominent feature in the top panels is the discrepant location of Class 2 (blue): In the right panel, the analytical solution, all blue circles are in the high ∆j z regime as expected (see section 3.4).However, in the left panel, the numerical solution -blue marks appear both in the high ∆j z regime and also mixed with Class 4 (cyan crosses) in the low ∆j z regime.This discrepancy arises due to the nature of the potential having a local maximum and the high sensitivity regarding whether E > V max or not.In the analytic solution, Class 2 represents a small region in the initial conditions phase space which is continuously connected to the region occupied by Class 4 but has a discontinuous outcome in δ min as mentioned.It is therefore unsurprising to have some misclassification between Class 2 and 4 leading to a notable deviation in ∆j z .As expected according to the arguments mentioned above the discrepancy appears in j max z and not in j min z as can be seen in the middle and bottom panels.
Finally, note that the tight correlation observed in the rightward tail in ∆j z is much tighter than in j max z and j min z separately, implying a correlation between them.While the analytical solution does not reproduce the tight correlation in ∆j z it is similarly narrower than the corresponding scatters in j max z and j min z in isolation.

SUMMARY
In this work, we have generalized the analytic solution developed in Klein & Katz (2023) for a precession rate β in the vicinity of β 0 .Using the analytic solution we obtain predictions for the amplitude of change in j z based on the initial conditions only, for librating KLCs.These analytic predictions agree approximately with the quantitative results of a full numerical solution of Equations 2-3.The analytic equations are equivalent to a particle moving in a one dimensional potential.We have shown that different features in the rich structures observed in the numerical results correspond to four classes of initial conditions: (i) No local maximum in the potential, (ii) Energy higher than the local maximum, and two classes with energy below the local maximum in one of two local minima (see section 3.4).
Strictly speaking, the effective equations (based on the ansatz in Equation 7) and analytic solution were developed close to the fixed point C K = −1.5 and for precession rates around β 0 (and for α ≪ 1 and j 0 z ≪ 1).In this regime, β 0 is the natural resonant frequency of the equations.Intriguingly, the generalized solution approximately agrees with the numerical outcomes for the majority of librating KLCs, even when C 0 K approaches values near 0, far from the fixed point and across a wide range of β values.
In future work we plan to analytically investigate the rotating KLCs (i.e where C K stays positive during the evolution).Numerically, the rotating class also shows rich structures and large changes in j z (not shown here).Given the success of the analysis around the fixed point for librating KLCs, we plan to try a similar approach for rotating KLCs by studying the vicinity of its fixed points at C K = 1 (specifically, j = 0, e z = 0, e x = cos (φ), e y = sin (φ) for any φ).Finally, there is a regime of initial conditions where C K changes its sign during the evolution.This class has low C 0 K and has numerically shown a lower amplitude of change in j z than the other two classes (for α ∼ 1 • ).For this class other methods should be used for developing an analytic understanding.

z
− j min z ) versus C 0 K is shown for various values of α and β in Figure 1 using black crosses.For each pair of α and β values the results are shown in one subplot.In the left column of panels, we maintain α = 1 • and explore different β values.In the right column of panels, we set β = 2.4 and examine different α values with equal logarithmic spacing.

Figure 1 .
Figure 1.Amplitude of the change in jz vs. initial CK (equal to the global constant in Equation 5 with j 0 z = 0 up to first order in α) for librating KLCs.Left panels: constant α = 1 • and different values of β (written explicitly above each subplot).Right panels: constant β = 2.4 and different values of α (logarithmically spaced).Results from a full numerical solution of Equations 2-3 are marked using black crosses.Green circles denote the prediction of the analytic solution using Equations 19, 22 and 23 (see section 3).The number N marked above each plot is the number of initial conditions shown.This number is the number of initial conditions that are librating KLCs among a random sample of 1000 realizations of initial conidtions.

Figure 2 .
Figure 2. The potential V (Equation 23) in blue and the constant energy E (Equation 22) in black for β = 2.4 showing the four distinct classes of (V, E) pairs: upper left panel: Class 1 -no Vmax, upper right panel: Class 2 -E > Vmax, lower left panel: Class 3 -E < Vmax and δ < δV max and lower right panel: Class 4 -E < Vmax and δ > δV max .

Figure 3 .
Figure3.Amplitude of the change in jz vs. initial CK in the numerical integrations (same as the black crosses in Figure1) color coded by the analytically specified nature of the pair (V, E) as defined in section 3.4: red -Class 1 (no Vmax), blue -Class 2 (E > Vmax), magenta -Class 3 (E < Vmax and δ < δV max ) and cyan -Class 4 (E < Vmax and δ > δV max ).Panels as in Figure1.

Figure 4 .
Figure 4.The maximal amplitude of change in jz vs. β for α = 1 • .Shown are results from full numerical solutions in librating KLCs (black crosses) and as predicted from initial conditions by the analytic solution (green circles).The results for β = 1.2, 1.8, 2.4, β0 ≈ 2.9 and 3.5 are the highest values attained by ∆jz in the corresponding panels of Figures 1.

Figure 5 .
Figure 5. Results for librating KLCs are shown for a case of α = 1 • and β = 2.4.Upper panels: Amplitude of the change in jz vs. initial CK .Middle panels: j max z vs. initial CK .Lower panels: j min z vs. initial CK .Left column: numerical results.Right column: analytic prediction.The color code in all panels is by the nature of the pair (V, E) as defined in section 3.4: red -Class 1 (no Vmax), blue -Class 2 (E > Vmax), magenta -Class 3 (E < Vmax and δ < δV max ) and cyan -Class 4 (E < Vmax and δ > δV max ).The top left panel is identical to the middle left panel of Figure 3.