Stability analysis of numerically exact time-periodic breathers in the Lugiato-Lefever equation: Discrete vs continuum

This paper reports an intimate connection between the recently observed frequency comb breathers in a continuous Lugiato-Lefever equation and breathers in an array of coupled resonators. Hopf and period-doubling instabilities and other features of the breathers associated to the Lugiato-Lefever equation are reported for the first time.


I. INTRODUCTION
Breathers are spatially localized nonlinear waves, with energy oscillating periodically in time. One of the paradigmatic continuous-wave equations for studying breathers in conservative systems is the integrable nonlinear Schrödinger (NLS) equation, with applications in many branches of physics (see, e.g., [1,2] and references therein). As these NLS breathers, known as Kuznetsov-Ma (KM) breathers [3,4], appear on top of a modulationally unstable nonzero background, they are unstable for infinite-size systems, but may be sufficiently stabilized for smaller systems, e.g., using periodic boundary conditions, so as to allow for their experimental observation in optical fibers [5].
Breathers also appear in time-periodically driven dissipative systems, typically arising from Hopf bifurcations of stationary dissipative solitons. Adding driving and damping terms to the NLS equation yields the Lugiato-Lefever equation (LLE) [6], which has become the paradigmatic equation for studying various nonlinear phenomena in optical cavities with dissipation [7][8][9]. In contrast to the conservative KM breathers, the nonzero background of the LLE breathers may be stabilized by dissipation, and there have been several experimental works in recent years confirming their existence as attractors for suitable parameter values and initial conditions, for macroscopic fiber cavities [10,11], and for optical microresonators [12][13][14][15][16][17][18].
On the other hand, there is currently also a large amount of interest in studying light propagation in spatially discrete systems [19], where the discretized version of the LLE appears as a model for coupled Kerr-nonlinear cavities [20][21][22]. Analogously to the cavity solitons of the continuous LLE, this model supports localized stationary modes also far away from its continuum limit, discrete cavity solitons, which were also shown to destabilize via Hopf bifurcations into oscillating solitons for certain parameter regimes [21].
It should be noted that previous numerical work investigating breathers in continuous as well as discrete LLEs was done by direct integration of the equations of motion, starting from suitable initial conditions yielding breathers as dynamical attractors. A disadvantage of this method is that it never converges to unstable breathers, and therefore it does not allow for continuation of breather families through regimes of instabilities. For example, breathers appearing from subcritical Hopf bifurcations of a stationary solution could not be followed to the bifurcation point. On the other hand, there are many existing continuation methods for calculating time-periodic orbits in dissipative systems to numerical precision by Newton-type iterations (see, e.g., [23]), and with a sufficiently accurate solution, linear stability analysis can be effectively performed through numerical Floquet analysis. This program has been carried out in many contexts for conservative systems and in particular used for determining the existence and stability properties of breathers in the discrete [24] as well as the continuous [1] NLS equation. For dissipative systems this type of analysis becomes numerically more challenging, as there are no longer continuous families of solutions with different time periods for given parameter values but rather solutions appearing at isolated points, with periods that have to be determined from the necessity to have balance between energy input and dissipation over one period. Nevertheless, a detailed investigation along these lines of regimes of existence and stability for time-periodic solutions of the parametrically driven and damped continuous NLS equation was reported in a series of papers [25][26][27]. However, here we perform the study for the LLE equation and explicitly for a lattice model, comparing solutions in the strongly discrete and continuum regimes.
In the present paper we employ a Newton algorithm for calculating time-periodic breathers of the discrete LLE to numerical precision and investigating their linear stability through Floquet analysis. Our aim here is not to provide detailed charts of existence and stability, as done for the parametrically driven damped NLS equation in [25,26], as this would require extensive computational resources. Rather, we wish to present the numerical framework and to show with a number of examples the relation between breathers in discrete and continuum systems. The outline of the paper is as follows. In Sec. II we present the continuous and discrete LLEs and discuss relations between their parameters. The numerical methods used for calculation and linear stability analysis of the breathers are described in Sec. III. In Sec. IV examples of our numerical results are presented: Sec. IV A shows the reproduction of the continuum breather stability regime for weak driving; Sec. IV B analyzes the continuation of the continuum breather towards weaker coupling; Sec. IV C shows the transition into the regime of period-doubling instability for the continuum breather with increasing driving; Sec. IV D illustrates the continuation of strongly discrete oscillating cavity solitons into continuum LLE breathers; and Sec. IV E confirms and explains the existence of stable LLE breathers in a regime of large detuning, which was recently observed experimentally [11]. Finally, some concluding remarks are presented in Sec. V.

II. MODELS
The dynamics of the temporal dissipative cavity Kerr solitons in microresonators can be described using the dampeddriven continuous nonlinear Schrödinger equation frequently referred to as the Lugiato-Lefever equation, with periodic boundary conditions: Here τ is rescaled time and the spatial variable φ measures the azimuthal angle of a ring resonator in a coordinate system moving with the group velocity [7,9]. Its discrete counterpart, modeling an array of coupled Kerr-nonlinear cavities, can be expressed as (with notation analogous to that in Ref. [21]) If viewed as a spatial discretization of (1) with N grid points spaced an angular distance φ, then n = φ φ , C = β 2 2( φ) 2 , and N = 2π φ . Thus, for sufficiently large C and sufficiently many lattice sites N, solutions to Eq. (2) (with γ = κ = 1) should approach solutions to Eq. (1) in a continuum limit, with → −α and E 0 → i f (the latter being equivalent to letting → iu; the phase of the driving term is arbitrary and can always be made real and positive by phase rotation of or u).
Note that Eq. (1) has been rescaled to have unity coefficients in front of dissipative and nonlinear terms (note though that the relative signs between these coefficients cannot be scaled away, but here have been assumed to be identical, from the particular physical setup considered). If the continuous equation (1) is considered for an infinite spatial interval, the dispersion coefficient β 2 may also be rescaled to unity by rescaling the spatial coordinate φ; however, for a finite system such a rescaling will change also the boundary conditions, and therefore we are left with three independent parameters β 2 , α, and f . Also, the discrete equation (2) may be rescaled to have κ = 1 and γ = ±1, with the free parameters C, , E 0 , and sgn(γ ). Explicitly, the rescaling amounts to setting C = C/κ, = /κ, E 0 = E 0 |γ |/κ 3 , u n = √ |γ |/κu n , and t = κt and dropping the primes. In addition, we may, by the staggering transformation u n → (−1) n u n , restrict consideration to C > 0.

III. NUMERICAL METHODS
There are many different, more or less sophisticated, ways to numerically obtain time-periodic solutions to continuous nonlinear dissipative wave equations such as (1). See in particular Ref. [23], which also contains many references to other earlier methods. Here, in order to directly connect results from continuous and discrete models, we take a slightly different route and instead use as a starting point the Newton method used previously to calculate breathing solutions in conservative and nondriven, discrete nonlinear Schrödinger equations. The method was briefly described in Sec. 4.1 of Ref. [24]. Here we need to amend the method since for the damped-driven system, there is no continuous family of timeperiodic solutions that can be continued versus the time period T , but rather isolated solutions with fixed, unknown time periods that have to be obtained as outcome from the Newton iteration. This feature was added with inspiration from works by Zemlyanaya and co-workers on the parametrically driven nonlinear Schrödinger equation (see, e.g., [25,26] and in particular [27] for details about their numerical scheme).
Briefly, we split the complex field u n into real and imaginary parts u n = x n + iy n and thus, for a system of N lattice sites, we obtain 2N coupled ordinary differential equations. We are looking for time-periodic solutions to (2) with unknown period T , u n (t + T ) = u n (T ). Since T is unknown, it is more convenient to rescale the time coordinate ast = t/T . We thus look for solutions with unity time period, x n (1) = x n (0) and y n (1) = y n (0), to the rescaled equations (below overdots indicate derivatives with respect tot) Linearizing the oscillation amplitudes around the exact timeperiodic solution x n → x n + ξ n and y n → y n + η n yieldṡ By simultaneous numerical integration from 0 to 1 of Eqs. (3) with initial conditions {x n (0), y n (0)} corresponding to a trial solution assumed to be close to the sought-after timeperiodic solution, i.e., fulfilling x n (1) ≈ x n (0) and y n (1) ≈ y n (0), and Eqs. (4) for 2N independent initial conditions for {ξ n (0), η n (0)}, we obtain the 2N × 2N Floquet matrix T 0 defined by The standard Newton iteration scheme [24] then amounts to solving the system of 2N linear equations for the corrections δu to the initial vector ( {x n (0)} {y n (0)} ) at each iteration step, and the norm of the vector on the right-hand side of (6) measures the accuracy of the obtained solution at that step. Now, to obtain also the additional unknown T , we consider variations T → T + δT in (3) and obtain, for the change of the amplitude vector at the final time [note that (3) is linear in T ], δx n (1) = ∂x n (1) ∂T δT = 1 Tẋ n (1)δT and δy n (1) = ∂y n (1) ∂T δT = 1 Tẏ n (1)δT . Thus, we add one additional column to the matrix 1 − T 0 in (6), and since our sought solution should fulfill alsoẋ n (1) =ẋ n (0) andẏ n (1) =ẏ n (0), the additional column can be taken as − 1 T ( {ẋ n (0)} {ẏ n (0)} ). We now have 2N equations for 2N + 1 unknowns. However, this is not a major problem since for any time-periodic trajectory there is always an additional degeneracy corresponding to the choice of origin of time. Thus, we may arbitrarily, according to our wishes, impose an additional condition, e.g., thatẋ n 0 (0) = 0 for some suitable chosen site n 0 or that x n 0 (0) = 0 (of course one has to be careful then that the chosen trajectory actually runs through the origin). Both these have been implemented and tested, and the (2N + 1) × (2N + 1) matrix can in general be inverted as required for the Newton iteration. However, it seems slightly preferable numerically to avoid adding this additional condition and instead use standard routines for singular-value decomposition to replace the matrix inverse with the pseudoinverse of the 2N × (2N + 1) matrix. (The particular scheme used here is a slight modification of a code originally developed in [28].) Having converged (typically after five to six iterations) to a solution with desired accuracy (typically 10 −11 with standard double-precision FORTRAN), linear stability analysis is then immediate by diagonalization of the matrix T 0 as obtained by integrating the linearized equations for the converged solution. Linear stability is equivalent to having all eigenvalues inside (or on, for marginal stability) the unit circle. For the numerical integration, we generally use a Bulirsch-Stoer algorithm, which empirically has been seen to perform well for these types of applications.
As we will mainly consider breather solutions that are spatially symmetric, we need only apply the Newton scheme to the N + 2 independent unknown variables, while for the Floquet stability analysis perturbations on 2N independent variables need to be considered. If the oscillating solutions are asymmetric in space (as some of the solutions discussed for the discrete model in [21]), the computational time will evidently be longer.

IV. RESULTS
For the continuous case, there are many papers showing various existence and bifurcation diagrams for solitons, breathers, etc., e.g., [29][30][31][32]. In particular, we may refer to Fig. 9 in [31], which indicates a regime B where timeperiodic localized breathers appear as stable attractors for the infinite LLE. A similar picture can be seen in Fig. 16 of Ref. [32] for the ring geometry. For the discrete case, the only paper that we are aware of explicitly showing existence and stability regimes for time-periodic solutions as attractors is [21], where Fig. 1 shows intervals in C for small and moderate values of C (not close to the continuous limit), only for one specific set of fixed values for and E 0 . Here we will show with a number of examples how the knowledge of the complete Floquet eigenspectrum helps to interpret and understand the dynamics in various regimes and how oscillating solutions in the discrete and continuous models may be connected.

A. Continuum regime: Continuation versus frequency detuning for weak driving
As our first example we show in Fig. 1 results from the stability analysis for a family of breathers close to the continuum limit for a rather weak driving, E 0 = 4, continued versus the frequency detuning. (In the continuum limit, the smallest possible value of the driving where localized breathers in an infinite system are observed is E 0 ≈ 2.65 [31].) In this regime, the breather family is always stable between its lower 033196-3 The lower branch corresponds to the family continued from larger C while the upper branches correspond to the bifurcating breather families, which are always unstable. Also shown are snapshots for C = 40 of (c) breather intensity |u n | 2 , and two corresponding eigenmodes of the linearized problem with eigenvalues on the negative real axis: (d) antisymmetric unstable and (e) symmetric stable (which will however destabilize for slightly smaller C). and higher limits of existence (at ≈ −8.5 and ≈ −4.8, respectively, for E 0 = 4) signaled by eigenvalue collisions at +1 in Fig. 1. The lower limit corresponds to the point where the breather Hopf bifurcates from the stationary cavity soliton, while the upper limit corresponds to a transition to a regime of spatiotemporal chaos as discussed in [31] (dotted line in Fig. 9 of Ref. [31]). The breather dynamics for these parameter values is analogous to Fig. 11 in [31].
Note that there is almost a period-doubling bifurcation around = −6, signaled by an eigenvalue approaching (but not reaching) −1; this is also consistent with Fig. 9 of Ref. [31], indicating that a period-doubled state would appear as an attractor for a slightly larger value of E 0 (cf. Sec. IV C). Note also the two eigenvalues that stay close to +1 during the continuation. One of them corresponds to the time translation along the periodic orbit and thus should exist for any time-periodic solution; its closeness to +1 can be used as a measure of the accuracy of the numerically obtained solution. Typically, its deviation from 1 is of the order of 10 −6 if the solution is obtained with accuracy 10 −11 . The other eigenvalue corresponds to the spatial translation, which is an exact continuous symmetry for the continuous equation but not for the discrete one. Thus, the closeness of this eigenvalue to +1 can be used as a measure of the accuracy of the discrete model in approximating the continuous one. Here, for C = 75 the maximum deviation is around 15%, which as shown below could be improved considerably by increasing C to approximately 100 (however, with the consequence that boundary effects increase due to the broadening of the solution with C, unless also the number of sites N is increased, which however increases the computational time considerably).

B. Continuation of continuum breather towards smaller coupling
As a second example, we show in Fig. 2 a continuation of the family of continuous breathers above for = −7 versus coupling constant C. First, as discussed above, note in Fig. 2(a) how the translational eigenvalue closely approaches +1 as C approaches approximately 100. Second, note how this localized eigenvalue, as C decreases, seemingly moves through the band corresponding to small-amplitude tail oscillations and reappears for C 50 on the negative real axis, where it passes through −1 and gives rise to a period-doubling bifurcation where a stable asymmetrically oscillating soliton should appear. This discreteness-induced instability seems analogous to the oscillations in the Peierls-Nabarro barrier found in [21] in a different parameter regime (see Sec. IV D below). Third, as C is further decreased an additional perioddoubling instability appears, this time corresponding to spatially symmetric oscillations as indicated by the eigenmode structure [ Fig. 2(e)]. Fourth, the continuation for decreasing C is lost at a saddle-node bifurcation approximately at C 21.96 [eigenvalue collision at +1 in Fig. 2(a)].
During the continuation towards smaller coupling, the time period slowly increases as shown in the lower branch of Fig. 2(b), where we also included a continuation of the bifurcating solution branches continued from the bifurcation point C 21.96. As can be seen, this snaking solution branch, which is always unstable, corresponds to breathers with larger period. It is also worth remarking that there is a tiny regime of stability for the lower branch very close to the bifurcation point, for 21.96 C 22, as seen from Fig. 2(a), where the two unstable eigenvalues again cross −1.

C. Increasing driving in the continuous regime
As a third example, we show in Fig. 3 the continuation for increasing driving strength E 0 of a breather in the continuouslike regime (C = 99) for = −7. As mentioned above, we now catch the period-doubling instability indicated in Fig. 9 of Ref. [31] roughly at E 0 4.5, and the corresponding oscillation mode is spatially symmetric around the breather center as shown in Fig. 3(c). As can be seen from the snapshots [ Fig. 3(b)], a further increase of E 0 results in a reshaping and considerable broadening of the (unstable) breather, which now more resembles a multibreather complex. Evidently, this state is more strongly affected by the boundaries. The continuation is finally lost in a bifurcation at E 0 = 6.2244 . . ., which again should be connected to the transition to the regime of spatiotemporal chaos [31], as discussed above for the continuation towards smaller frequency detuning.

D. Continuation of strongly discrete breathers
As a fourth example we show in Fig. 4 the continuation versus coupling constant C for a highly discrete breathing soliton, with parameter values chosen as in Ref. [21]. This breather Hopf bifurcates from the stationary discrete soliton at C 2.3, as shown in Fig. 1 of Ref. [21]. The period-doubling instability here appears around C 7, which coincides with the point denoted by C B3 in Ref. [21]. As can be seen in Fig. 4(c), the corresponding unstable eigenmode is antisymmetric around the breather center, which causes its motion in the lattice as found in [21]. As indicated in [21], this is related to the period-doubled solution itself being unstable at this parameter value. The breather restabilizes again at C 10.5 (point C B1 in [21]), slightly before it again merges with the stationary discrete soliton at C 10.7 [eigenvalue reaching +1 in Fig. 4(a)].
Note that the parameter values used above for the strongly discrete system, E 0 = 2.0 and = −3.0, correspond to driving and frequency detunings below the lower limit for the existence of breathers in the continuum regime, which according to Fig. 9 of Ref. [31] are E 0 ≈ 2.65 and | | ≈ 3.5. To search for a strongly discrete breather which can be directly continued to a continuum breather, we first continue the solution above towards larger E 0 and | |. As shown in Fig. 5, for E 0 = 3.0 and = −4.0 there is indeed a smooth continuation in C from C ≈ 2.4, where the breather Hopf bifurcates from a strongly discrete cavity soliton and then smoothly transforms to a continuumlike soliton for large C, indicated by the translational eigenvalue smoothly approaching 1 for C 50. Comparing Fig. 5 with the corresponding continuation of the continuumlike breather for larger driving and detuning in Fig. 2, the instability scenario is qualitatively similar with two period-doubling instabilities corresponding to symmetry-breaking and symmetry-conserving eigenmodes, respectively. However, while the continuation towards smaller C in Fig. 2 is interrupted at the bifurcation point at C ≈ 21.96 where the slope of the curve showing the time period versus C becomes infinite, the corresponding curve for the smaller driving and detuning in Fig. 5(c) instead reaches a maximum T ≈ 2.92 for C ≈ 10.6 and then again decreases until the strongly discrete breather reaches the point where it bifurcates from the stationary discrete cavity soliton.
Thus, to summarize this section, a major conclusion is that the oscillating discrete cavity solitons found in [21] are indeed lattice versions of the breather solitons of the continuous LLE, as there is a smooth continuation in parameter space connecting them.

E. Breathers at large detunings
Our last example considers a regime of large detuning and driving where, somewhat unexpectedly, stable breathers were found recently, experimentally as well as in direct numerical simulations for the continuous LLE [11]. With our notation, the regime studied in [11] corresponds to E 0 = 20 √ 5 ≈ 45 and −90 −80. However, as can be inferred, e.g., from Fig. 13 of Ref. [29], for such large values of E 0 , stable stationary cavity solitons should exist for −π 2 E 2 0 /8 −(E 0 /0.08 √ 2) 2/3 , where the lower limit corresponds to the existence threshold [33] and the upper limit to the instability threshold through Hopf bifurcation. Thus, for E 0 = 45, stable stationary solitons exist for −2500 −54 and so there is a regime of coexistence between the stable stationary soliton and stable breather, the origin of which was unclear in Ref. [11].
Attempting to reproduce the breather from [11] within the framework of the discrete model (2), we first note that very large values of C are needed in this regime to get close to a continuumlike breather, due to its large amplitude. Specifically, with = −85 we could obtain a breather for C = 500, but suffering from discreteness-induced period-doubling instabilities analogous to those discussed in Secs. IV B and IV D. Increasing to C = 700 a stable breather is indeed obtained, and below we will illustrate the continuation of this breather versus detuning. However, the consequence of using such large values of C is that the solution will broaden significantly, and finite-size effects will no longer be negligible.
Thus, we stress that, although we obtain for our rather small finite system a similar existence region for a stable breather as in [11], the detailed nature of the bifurcations limiting its existence may differ for larger systems.
In Fig. 6 we illustrate the continuation of the breather found for E 0 = 45, C = 700, and N = 121 sites versus frequency detuning. For these parameter values, the existence regime is bounded by two saddle-node bifurcations, at ≈ −96.93 and ≈ −54.75. The breather remains stable for −81.6, where a weak Hopf instability appears, and it becomes strongly unstable through period doubling when −80.9, as illustrated by the large negative (real) eigenvalue in Fig. 6(a). The solutions bifurcating with the breather at the saddle-node bifurcations [blue and green lines in Fig. 6(b)] are always unstable. As can be seen from Fig. 6(b) (purple line), the time period for the stable breather increases monotonically with increasing and continues to increase also in its unstable regime. In the inset of Fig. 6(b) we have superposed the plot of the time period around the upper bifurcation point with a plot of twice the period of the stable fundamental breather, which was found to exist for −56 −52 for the parameter values used. Thus, although the solution continued from larger negative apparently does not directly bifurcate with the stable breather, the relation between their periods approaches 2 in this regime.
As a final illustration, we show in Fig. 7 the resulting dynamics for a few oscillations of breathers at different locations in Fig. 6. Figure 7(a) shows the stable breather at large negative , which essentially reproduces the dynamics from Fig. 4 in Ref. [11]. Notice that each amplitude peak is in fact a double peak. Moving into the unstable regime of the purple branch in Fig. 6(b), the period increases and the double-peaked feature is lost, as shown in Fig. 7(b).   Fig. 6(b)] with a period very close to half that of the breather on the purple branch. As can be seen, this solution has much smaller oscillation amplitude and much smaller peak amplitude and, in contrast to the other solutions which have visible oscillations also in the tails, its oscillations outside the central region are negligible. Finally, Fig. 7(d) corresponds to the unstable solution on the blue branch of Fig. 6, which is characterized by considerably larger tail oscillations. A similar scenario was also seen for the parametrically driven and damped NLS equation (see Fig. 3 of [25]). Note that the regime of large detuning and driving studied here is, by the rescalings mentioned in Sec. II, equivalent to a regime of small damping. Thus, one may anticipate that some of these solutions could be considered as dissipative counterparts to breather solutions of the conservative NLS equation. In [16], the KM breather of the NLS equation was used to derive an approximate ansatz for a dissipative LLE breather, for the case when the amplitude oscillations of the background are small compared to the main breather peak amplitude. A distinguishing feature of this ansatz is that the dissipative breather frequency becomes identical to the detuning, i.e., T = 2π/| |, which was also seen in [16] to be close to (but slightly larger than) numerically as well as experimentally determined frequencies. For our solutions illustrated in Figs. 6 and 7 we observe the same tendency for the time period of the fundamental stable breather [ Fig. 7(c)], i.e., T 2π/| |, while the periods of the other solutions are approximately twice this value. Whether there indeed exists a smooth continuation in parameter space between the conservative KM breather and the fundamental dissipative LLE breathers is a more challenging numerical issue, which is left for future work.

V. CONCLUSION
We have implemented a Newton-type algorithm for calculating time-periodic breather solutions to the discrete LLE and used numerical Floquet analysis to investigate the linear stability properties of the breather families resulting from continuation versus frequency detuning, driving, and coupling constant, in various regimes of parameter space. A major aim of this work has been to unify previously existing results for the discrete and continuous LLEs into a common framework, and as we find a smooth continuation between strongly discrete and continuumlike breathers with increasing coupling constant, we may conclude that the oscillating discrete cavity solitons from [21] are indeed lattice versions of the continuum LLE breathers. In general, having access to the full spectrum of Floquet eigenvalues provides a deeper understanding of the various instability mechanisms and bifurcations than earlier studies relying on direct numerical simulations could provide. As a particular example of this, we considered the regime of large detuning and driving with coexistence between stationary solitons and breathers recently found experimentally in [11]. Our results (for a finite ring with large coupling constant) not only confirm the existence of linearly stable breathers in this regime, but also show the fate of this family when continued versus detuning: Its existence regime is limited by two saddle-node bifurcations with unstable breathers and its stability regime is limited from one side by Hopf and period-doubling instabilities.
As pointed out in the Introduction, obtaining a detailed existence and stability chart for breathers in the full threeparameter space of the discrete LLE would require extensive computational efforts and is left for future work. We have also left a number of other interesting issues unexplored. For example, it would be straightforward to use the same techniques to analyze also the properties of asymmetrically oscillating breathers; however, as we pointed out, computational times would increase considerably. Another interesting issue would be to systematically follow the sequence of period-doubled oscillating breathers with increasing periods; however, this also would require longer computational times. Moreover, an important issue is the dependence of system size. Here we chose for numerical convenience a rather small system with 121 sites; this could certainly be enlarged but again with the consequence of rapidly increasing computational times. In the strongly discrete limits, we are confident that our results approximate those of an infinite system with good accuracy, but as we indicated, there are several regimes close to the continuous limit where the oscillations near the boundaries are considerable. Obtaining reliable predictions for the infinite systems in these regimes would require studying some sequence of increasing lattice sizes; however, these issues are left for future consideration.