HarmonicBalance.jl: A Julia suite for nonlinear dynamics using harmonic balance

HarmonicBalance.jl is a publicly available Julia package designed to simplify and solve systems of periodic time-dependent nonlinear ordinary differential equations. Time dependence of the system parameters is treated with the harmonic balance method, which approximates the system's behaviour as a set of harmonic terms with slowly-varying amplitudes. Under this approximation, the set of all possible steady-state responses follows from the solution of a polynomial system. In HarmonicBalance.jl, we combine harmonic balance with contemporary implementations of symbolic algebra and the homotopy continuation method to numerically determine all steady-state solutions and their associated fluctuation dynamics. For the exploration of involved steady-state topologies, we provide a simple graphical user interface, allowing for arbitrary solution observables and phase diagrams. HarmonicBalance.jl is a free software available at https://github.com/NonlinearOscillations/HarmonicBalance.jl.

Despite the ubiquitous presence of nonlinear dynamical systems in science and engineering [ Fig. 1b], their behaviour is not analytically tractable in most cases. Hence, one often resorts to using numerical ODE solvers e.g. as in molecular dynamics simulations [17]. These usually focus on initial value problems, where the system's state is advanced from a set of initial conditions. Linear systems, given sufficient time to freely evolve, usually relax to a unique stationary or steady state, i.e., a state where the system no longer evolves in time. In a harmonically driven system, the steady states will typically display time dynamics which are also harmonic, requiring a corresponding definition of a dynamical stationarity condition.
Nonlinear systems challenge this approach since they can feature multiple steady-state solutions. The knowledge of all such solutions is of key importance for analysing the system's behaviour, revealing experimentally observable phenomena such as hysteresis, spontaneous symmetry-breaking, or noise-induced switching dynamics [18][19][20]. However, when using a dynamical solver (or a time-resolved experiment), only a single steady state is found per run, depending on the initial conditions [21,22]. Therefore, a complete exploration of the solution landscape would require infinite sampling of the continuous space of initial conditions. Basic functionality of HarmonicBalance.jl: A nonlinear (polynomial) ODE system with harmonic time-dependence is defined along with a set of expected response frequencies for each resonator, ω i j (red input boxes). The package finds all dynamical steady states, identifies the most likely response given a parametric sweep, and supports numerical time-evolution within the harmonic ansatz (introduced in 2.1).
An alternative approach to finding steady-state solutions in harmonic systems is transforming the system into a frame rotating with the applied drives [23,24]. In such a frame, the drives and corresponding steady states appear stationary, reducing the problem to finding the roots of a time-independent system of nonlinear equations [25,26]. While this approach is standard in few-variable problems, the proliferation of roots in multivariate nonlinear systems constitutes a challenge in itself [27][28][29]. Indeed, the most straightforward root-finding algorithms (e.g., the Newton-Raphson method) only find a single solution in the vicinity of an initial guess, making the exploration of all steady states infeasible.
Fortunately, many of the aforementioned nonlinear examples [cf. Fig. 1b] feature ODE systems with polynomial dependence on the dynamical variables and their derivatives. For these, seeking all steady states in a given rotating frame generates coupled polynomial (algebraic) equations, which are numerically solvable using the homotopy continuation method [30][31][32]. This method utilises the continuous deformation of an exactly solvable problem into the problem of interest, thus finding all the roots of coupled polynomials in one go. To our knowledge, however, no homotopy-continuation-based toolbox for the analysis of steady states of harmonically-driven systems has been developed so far.

Harmonically-driven nonlinear systems: basic principles
HarmonicBalance.jl focuses on harmonically-driven nonlinear systems, i.e., dynamical systems governed by equations of motion where all explicitly time-dependent terms are harmonic. In the exposition here, we will assume a general nonlinear system of N second-order ODEs 1 with real variables x i (t), with index i = 1, 2, · · · , N and time t as the independent variable, The vector x(t) = (x 1 (t), ..., x N (t)) T fully describes the state of the system. Physically, x(t) encompasses the amplitudes of either point-like or collective oscillators (e.g., mechanical resonators, voltage oscillations in RLC circuits, an oscillating electrical dipole moment, or standing modes of an optical cavity). We assume F(x(t), t) can be decomposed into a sum of L periodic terms 2 :, with vector fields f l (x). The field f 0 describes static properties of the system while f l =0 account for explicit time-dependence, typically induced by one or more sources of periodic driving and/or parameter modulation with frequencies {ω l } and phases {φ l }. In Table 1, we list several examples of terms commonly constituting Eq. (2). Note that in a linear system, F(x, t) can be written as F(x, t) = M(t)x + b(t), where the matrix M(t) contains spring constants and linear couplings, while b(t) is a vector of external forces. Diagonalising M(t) then yields the so-called normal modes of the system. While the notion of normal modes does not directly apply in a nonlinear system, they usually constitute a convenient basis choice for a perturbative treatment.

The harmonic expansion
For sufficiently long times (i.e., after any transient responses have disappeared), the solutions of Eq. (1) are expected to appear as a sum over harmonics. Let us illustrate this point using the example of driven harmonic oscillators. In this simple case, Eqs. (1) and (2) take the form where we assumed a constant vector g and drive frequency ω d . The standard method to solve for the steady states is to Fourier-transform both sides of Eq. (3); the resulting equations give an exact solution for x(t) in terms of its Fourier coefficients, (2) and their origins. The rightmost column shows the frequency conversion taking place, assuming two variables x i , x j , oscillating at frequencies ω i , ω j (for example, ω i → 2ω i means, that oscillating at frequency ω i leads to additional oscillations at frequency 2ω i ).
is the Dirac delta function of z centred at z 0 . This procedure is effective because in the Fourier domain, the l.h.s. of Eq. (3) becomes diagonal (i.e., it involves a single frequency) while the applied drive reduces to a Dirac delta function centred at ω = ω d . We can thus select a single frequency out of the continuous space of all frequencies for which the solution is nonvanishing. Furthermore, due to the linear superposition principle, responses to arbitrary driving terms can be constructed out of the solution in Eq. (4).
The same task becomes intractable if nonlinear terms are introduced. Nonlinearities facilitate frequency conversion by coupling different harmonics of the system, rendering F(x, t) non-diagonal in Fourier space. As an example, let us consider a single driven Duffing oscillator [77], whose equation of motion reads where ω 0 is the natural frequency, α is the nonlinear coefficient, F is the drive amplitude, and x is now a scalar. The nonlinear (Duffing) term in Fourier space reads coupling all combinations of four harmonics that sum to zero. This results in frequencyconversion processes known as four-wave mixing, which here convert, to lowest order, the driven oscillation at frequency ω d to frequency 3ω d ; four-wave mixing appears as off-diagonal terms in Fourier space. The frequency conversion subsequently propagates through the entire spectrum, generating an infinite number of Fourier components. A nonlinearity thus precludes a closed-form solution of the problem in Fourier space; this is a common trait of nonlinear physical systems; see Table 1 for examples of frequency-converting effects.
A serviceable approach to numerically approximate the harmonics of a driven nonlinear system involves truncating the spectrumx(ω) to a finite set of frequencies [78]. The idea of truncation in Fourier space is at the core of several widely-used methods, such as harmonic balance [79,80], the rotating-wave approximation [81], the van der Pol transformation [4] in combination with Krylov-Bogoliubov averaging [82], Magnus expansion [83,84], secular perturbation theory [77], and also appears in the contemporary concept of Floquet engineering [85,86]. We implement the approach of harmonic balance using a generalised ansatz [87,88] where the sum runs over the finite set of desired frequencies {ω i, j } describing the coordinate x i . Here, T represents a "coarse-grained" timescale that is much slower than the oscillations in the system (T 2π/ min{ω j }). Equation (7) represents an attempt to capture the dynamics of the system using a discrete set of real functions {u i, j (T ), v i, j (T )}, which in our package we dub the harmonic variables, while the fast oscillations are accounted for by the sine/cosine terms. Note that this ansatz would be exact if the set of frequencies {ω i, j } covered all real numbers. Generally, it is not straightforward to identify the relevant frequencies {ω i, j }. A good starting point is the frequencies of any external drives, combined with the frequency conversions presented in Table 1. Alternatively, one can pick the highest-weight discrete Fourier components in a time trace obtained numerically or from experimental data.
Plugging the ansatz Eq. (7) into the ODE (1), we obtain equations governing pairs of harmonic variables (u i, j , v i, j ) by convolving both sides with cos ω i, j t and sin ω i, j t , respectively (see Appendix A). The harmonic variables themselves are treated as constants during this step. We thus obtain two equations for each ω i, j . To reflect the slowly-changing nature of the harmonic variables, we drop all of their time derivatives of order two or higher, e.g., u i, j ,v i, j → 0. The resulting set of equations has no explicit time dependence; we call these the harmonic equations,  [18,82] and the rotating-wave-approximation [81]. Note that our scheme scales up very fast -each harmonic ω i, j of a variable x i is converted into two harmonic variables. Hence, for a system of N interacting components, each expanded in M harmonics, Eq. (8) consists of 2N M harmonic equations. Note that the discussed approach is not exclusive to fixed points, and the ansatz can account for periodic asymptotic orbits, i.e., to limit cycles [1-5, 41, 89-91]. Limit cycles can emerge when the system loses stability around Hopf bifurcations and a selfsustained oscillation emerges. We plan to include limit cycle functionality in future releases of the package.

Solving the harmonic equations for steady-state solutions
There are, in principle, two ways to extract the long-time behaviour of the system from Eq. (8). First, one can make use of an initial-value ODE solver to numerically propagate a state in T , starting at u(T = T 0 ). Given sufficient time, the system will typically converge towards a steady state u 0 . Secondly, we can look for steady states by explicitly requiring the l.h.s. of Eq. (8) to vanish, leaving us with the task of solving the nonlinear algebraic equations system The second approach is advantageous since, in principle, all steady states of the system can be found, irrespective of whether or how they can be reached by time evolution. In general, however, obtaining the roots of coupled nonlinear algebraic equations is highly non-trivial: for n coupled polynomial equations of order p, Bézout's theorem [92] provides an upper bound p n R e ( s )°2°1 to the roots of the harmonic balance equations [cf. Eq. (9), Appendix A] (green dots), denoted u 0 = (u 0 , v 0 ), using homotopy continuation. In the algorithm to track multiple parameter steady states, we carry out a two-step process encompassing a warm-up tracking (deformation parameter λ w ∈ (0, 1)), where singular paths (not-shown) are filtered-out, followed by a subsequent tracking to the harmonic balance equations, (deformation parameter λ p ∈ (0, 1)) [see Appendix B for further details]. The roots u 0 traverse the complex plane until paths intersect with the planes Im(u 0 ) = 0 and Im(v 0 ) = 0 (light gray), at the 3 real roots of F(u 0 ). For this figure α = 1, ω d = 1.03, ω 0 = 1, F = 0.01, θ = 0, and γ = 0.01.
for the number of distinct complex solutions. For example, in the Duffing oscillator shown in Eq. (5), we have one variable x -given an ansatz Eq. (7) with a single harmonic. This generates 2 polynomial harmonic equations of order 3, leading to a maximum of 3 2 = 9 solutions. Generally, for N coupled Duffing oscillators each expanded in M harmonics, we have 2N M harmonic equations of order 3, leading to an upper bound of 3 2N M solutions. The exponential scaling rapidly makes charting the complete steady-state solution landscape a formidable challenge. Accordingly, many applications of harmonic balance deal with few-variable systems and rely on the perturbative treatment of nonlinearities.
In our package, we solve the algebraic Eq. (9) using the homotopy continuation method [93,94] as implemented by the open-source package HomotopyContinuation.jl [31]. To find the roots of a polynomial, this method, in its simplest form, starts from another analytically-solvable polynomial of the same order, which saturates the Bézout bound. For instance, to find the roots of a single, p-th order polynomial P(z) of the variable z, one can start from the polynomial U(z) = z p − 1 with roots e 2πik/p (k = 1, ..., p). Employing this so-called total-degree homotopy, the polynomial U(z) is "slowly" deformed 3 into the polynomial P(z), correcting the known roots (i.e., tracking the roots) after each step. At the end of the homotopy continuation, the obtained set of roots is guaranteed to be complete, as it has been tracked from the full set of p complex roots of U(z).
A generalisation of the total degree homotopy approach enables the solution of Eq. (9), by instead tracking the roots of the uncoupled system U(u 0 ) = {u 2 −1, · · · }, with d r equal to the degree of the polynomialF r (u 0 ), to those ofF(u 0 ) (here r = 1, 2, · · · , 2N M is an index labelling all harmonic balance equations). Such an approach leads to a number of solution paths to track that scales exponentially with N M . This constitutes a challenge in usual exploratory research scenarios that seek to find the behaviour of the solutions as parameters vary. Crucially, in physical problems (and in contrast with single polynomials),F(u 0 ) often displays a vastly smaller number of roots than the Bézout bound (see Examples). This implies that many solution paths from U(u 0 ) = 0 are usually irrelevant (singular 4 ) since they do not lead to roots ofF(u 0 ). In HarmonicBalance.jl, we harness a common strategy to reduce the computational overhead in multi-parameter root finding [30]. This strategy consists of the two-step algorithm described in Appendix B and exemplified in Fig. 2.
Finally, note that although the obtained solutions are complex, only real roots of Eq.(9) are physically meaningful 5 . For a more detailed description of the method and its implementation, see the documentation of HomotopyContinuation.jl and references therein [30,31].

Stability analysis and linear response
Let us assume that we found a real solution u 0 of Eq. (8). When the system is in this state, it responds to small perturbations either by returning to u 0 over some characteristic timescale (stable state) or by evolving away from u 0 (unstable state). To analyze the stability of u 0 , we linearize Eq. (8) around u 0 for a small perturbation δu = u − u 0 to obtain where J(u 0 ) = ∇ uF | u=u 0 is the Jacobian matrix of the system evaluated at u = u 0 6 .
The solution to Eq. (10) can be expanded in terms of the complex eigenvalues λ r and eigenvectors v r of J(u 0 ), namely The dynamical behaviour near the steady states is thus governed by e λ r T : if Re(λ r ) < 0 for all r, the system returns to u 0 under a small arbitrary perturbation and the state is classified as stable. Conversely, if Re(λ r ) > 0 for any r, the system is unstable -perturbations such as noise or a small applied drive will force it to evolve away from u 0 .
As the system is coupled to dissipative baths, it will be also subject to fluctuations. While the current focus of our package is deterministic dynamics, analysis of weak fluctuations around a stable fixed point is included. The linear response of a stable, steady state to an additional oscillatory force, caused by weak probes or noise, is often observed in experiments [95]. It can be calculated by adding a small driving term δf cos(Ω d T ) to the harmonic Eq. (8). To account for the time dependence of the perturbation, linearisation around u 0 should retain the previously-dropped higher-order time derivatives. While implemented in the package, we leave a thorough discussion of this topic, as well as noise-activated dynamics, to future work.

Structure of HarmonicBalance.jl
The bulk of HarmonicBalance.jl is written in Julia, a language combining the accessibility of interpreted languages such as Python with the performance of compiled languages such as C and Fortran. Our package is designed with simplicity of use and scalability in mind. We give a short overview of the package structure and basic workflow below; it is also displayed in Fig. 3.
1. ODE systems serve as primary input and are inserted in symbolic form via Symbolics.jl.
Non-autonomous harmonic systems are accompanied by a user-defined set of frequencies to build the ansatz [cf., Eq. (7) and Fig. 1c]. From this, a set of harmonic (autonomous) equations is obtained symbolically.
3. Steady-state solutions are further analysed, e.g., for stability and other criteria. A flexible toolbox for the calculation of observable quantities is included, with a user-defined mapping of steady-state solutions and the calculation of linear response spectra.
4. Steady-state solutions for one-and two-dimensional parameter sets are readily visualised. Interactive plotting routines help exploratory analysis of the rich topology of the solutions (e.g., bifurcations) for multi-dimensional parameter sets. An ODE solver may be used to obtain the slow-time (T ) dynamics and verify the steady-state results.
In the following, we detail the working principles of HarmonicBalance.jl. For detailed online documentation on the latest release, we refer the reader to [96].

Defining a system, extracting harmonic equations.
Conceptually, specifying a system requires two ingredients: its ODE of motion, such as Eq. (1), and the set of oscillating frequencies forming the harmonic ansatz in Eq. (7). Once defined, the equation of motion is stored in the dedicated object DifferentialEquation. For this primary input and subsequent symbolic manipulations, we employ the Julia package Symbolics.jl [97], whose emphasis on high performance is essential in dealing with complex problems, such as that shown in Section 5.4.
After constructing the harmonic ansatz Eq. (7) and introducing the slow time T , the harmonic equations governing the harmonic variables u are found. We assume slow or no evolution in T and therefore drop any second and higher-order derivatives with respect to T as well as products of du d T . The equations are then Fourier-transformed, returning the harmonic Eq. (8) (for a concrete example, see Appendix A and Section 5.1). These equations are stored in the object HarmonicEquation, which in itself stores instances of HarmonicVariable, each specifying one pair {u i, j , v i, j }.

Obtaining and characterising steady states
To find and analyse the steady states of a HarmonicEquation, we need the corresponding (algebraic) steady-state equations and the Jacobian matrix; these are symbolically stored in the object Problem.
As the next step, numerical parameter values are specified. This is the required input for our primary solving method, get_steady_states, which allows the user to specify which parameters are constant and which are varied. One may vary any number of parameters, i.e., solution sets of any dimension are supported. We then employ a homotopy to retrieve steady states for each parameter set value. The native multi-threading support of HomotopyContinuation.jl can dispatch path tracking over multiple cores. Note that complex roots are being followed throughout the procedure -real solutions are only filtered out a posteriori.
Once steady states are found, HarmonicBalance.jl automatically classifies each state by whether it is real (complex roots bear no physical meaning) and by its stability. The final output is a Result object.

Visualisation
Functionality for static plotting of steady-state and time-dependent solutions is provided using the Matplotlib library [98]. At the core of our plotting routines is the function transform_solutions, which evaluates a symbolic expression by substituting each of the solutions from a Result object. This enables 1D and 2D plots of each solution and functions derived from them. Functionality for selection of transformed solutions (e.g. those with the maximal amplitudes) is available through the method map_multi_solutions 7 . We implemented classification functionality to help navigate potentially involved solution landscapes to show/hide/label solutions satisfying a given condition.
Rather than the solutions themselves, a phase diagram is often desired, which shows how the qualitative behaviour of the system changes in parameter space. We provide the functionality to distinguish parameter space regions by the number of (stable/all) solutions. We include an additional test that is able to detect qualitative differences between points in parameter space where the number of stable/unstable solutions is invariant (e.g. transcritical bifurcations) 8 .
Finally, we highlight the function interactive_phase_diagram_2D; this produces an interactive window where the solution landscape can be explored along a given dimension by simply clicking on different parts of a two-dimensional diagram.

Time-dependent simulations
The behaviour of a system can be found using a numerical ODE solver and an initial condition u(T = T 0 ). This approach is beneficial for analysing systems whose parameters are being adiabatically varied in time, i.e. slower than characteristic system reaction times -a common experimental approach to explore the solution landscape, detecting bifurcations and hysteresis. Other potential uses of time-dependent simulations are verifying steady-state results, their stability and fluctuation spectra, and identifying their basins of attraction.
A time-dependent simulation may use either a DifferentialEquation or a HarmonicEquation. The former represents the system in the time domain and uses no 7 In general, this function provides support for many-to-one functions of multiple transformed solutions. 8 In short, this function first labels each point in parameter space with a bit string indicating the stability of each solution there (e.g. [is_stable(sol_1),is_stable(sol_2),· · · ]=[011· · · ]). Next, each bit string is assigned to a unique integer identifier.
approximations. However, the numerical propagation of oscillatory dynamics can be extremely inefficient. Specifically, time grids need to resolve fractions of the period of the fastest oscillation to keep numerical precision. A HarmonicEquation is significantly faster to solve since it describes the system using the slow time T , where the oscillatory steady states appear time-independent. This approach, however, can only capture the chosen set of harmonics and nearby frequencies through the slowly varying amplitudes u. To include additional frequencies, the harmonic ansatz must be expanded. To illustrate the interface to time-dependent solvers, we provide a simple example in the online documentation https: //nonlinearoscillations.github.io/HarmonicBalance.jl/stable/examples/time_dependent/.

Comparison with other harmonic balance implementations
Steady-state problems in nonlinear periodic ODEs appear in diverse areas of science and technology. Several free and commercial packages exist to solve them. Examples of opensource software include Xyce -a high-performance parallel electronic simulator that can perform harmonic balance analysis [99]. The harmonic balance method is also natively supported for general nonlinear multiphysics finite element simulations in the open-source C++ FEM library Sparselizard [100]. An example of a commercial software package is Cadence AWR Microwave Office, which focuses on the analysis of electrical circuits in the frequency domain (i.e. calculation of voltage and current for a given set of elements), and includes features that are not yet implemented in HarmonicBalance.jl (e.g., noise analysis). Generally, these packages are use-case specialised, making their application in other settings challenging.
Likely the closest existing tool to our Julia package is NLvib [79]. However, it focuses on predefined problems primarily relevant to mechanical engineering and requires a MATLAB license to use, which is unaffordable for many research groups and educational institutes. Crucially, in all of the above packages, the harmonic balance equations are solved by either i) time evolution from a given set of initial conditions or ii) single-root finding methods such as Newton's. These approaches only return one steady state at a time, and even when combined with the continuation of a known solution, disconnected solution branches remain hidden. Our distinction from existing work lies in three main points: • The use of homotopy continuation allows us to find all possible solutions of the harmonic balance equations, where solutions for multiple parameters can be obtained reliably and efficiently. In this critical step, usually the most computationally intensive, we rely on HomotopyContinuation.jl, a package which outperforms other established homotopy continuation libraries [31].
• Arbitrary equations of motion can be entered and processed. This is crucial for use by the academic community, where an ab-initio approach to physical problems is often preferred over specialised GUI-based tools.
• The code and its dependencies are entirely open-source. This enables a natural synergy with Julia's existing rich ecosystem for scientific computing.

Examples of HarmonicBalance.jl usage
In this Section, we present several example problems. Instructions to install HarmonicBalance.jl and further detailed examples can be found at https://github.com/NonlinearOscillations/ HarmonicBalance.jl and in the links therein.  The simplest use case for HarmonicBalance.jl is a driven Duffing resonator governed by Eq. (5). For a driving frequency ω d in the vicinity of the natural resonance frequency ω 0 , the dominant behaviour can be captured by a single harmonic, the excitation frequency ω d . The basic lines of code follow below: definition of the system, implementation of the harmonic ansatz

Duffing oscillator
and derivation of corresponding harmonic equations proceed as using HarmonicBalance @variables α , ωd , ω0 , F , t , γ , x ( t ) # declare constant variables and a function x ( t ) diff_eq = D i f f e r e n t i a l E q u ation ( d ( d (x , t ) ,t ) + ω0^2 * x + γ * d (x , t ) + α* x^3~F * cos (ωd * t ) , x ) add_harmonic !( diff_eq , x , ωd ) # specify the ansatz x = u ( T ) cos The steady-state amplitude corresponding to the harmonic ω d is given by X ω d = u 2 1 + v 1 ; this can be inspected via plot ( solutions , x = "ωd " , y = " sqrt ( u1^2 + v1^2) " ) which produces curves similar to Figs. 4a,b (top row). In such curves, markers are employed to denote the stability of solutions, determined from the Jacobian matrix eigenvalues, which are evaluated at each point [ Fig. 4c]. Real and imaginary parts of Jacobian eigenvalues are displayed via a call to p l o t _ 1 D _ j a c o b i a n _ e i g e n v a l u e s ( soln_1d , x = "ωd " ) ; Moreover, simultaneous parametric sweeps over many parameter dimensions can be produced by adding ranges to swept, e.g., to get a 2D dataset, swept = ( F = > LinRange (0.0001 ,0.0041 ,150) , ωd = > LinRange (0.97 , 1.03 ,150) ) A 2D phase diagram depicting the total number of solutions, irrespective of stability [ Fig. 4d], can then be produced by plo t_2D _ p h a s e _ d i a g r a m ( solutions , stable = false , observable = " nsols " ) The cubic nonlinearity leads to multiple frequency upconversion processes (Section 2.1). In particular, we focus on the conversion of a resonantly excited oscillation at ω d ≈ ω 0 to an off-resonant oscillation at 3ω d 9 . The analysis of this mechanism is helpful to study convergence studies of the ansatz. To corroborate it, we implement the ansatz including harmonics at ω d and 3ω d : In Fig. 4 bottom, the nonlinearity, yields a finite amplitude X 3ω d = u 2 2 + v 2 2 > 0 for the third harmonic 3ω d = ω 0 . Comparison of phase diagrams for the two ansatzes in Fig. 4d does not reveal new regions in terms of the number of solutions. As described in Appendix A, this insight justifies a perturbative treatment of the upconverted response.

Navigating involved solution landscapes: Coupled Duffing resonators
The power of the homotopy continuation method lies in its ability to capture all possible solutions as parameters are varied. In contrast, experiments in nonlinear systems often implement parameter variations by adiabatic parameter sweeps, which can only access a subset of steady states. This is because such protocol reaches only one steady state at a time.
The knowledge of all possible solutions can be beneficial to interpret experiments of coupled nonlinear oscillators. As the number of oscillators increases, so does the dimensionality of the problem and the number of potential steady states (Section 2.2) with corresponding intricate solution topology. As an example of how increasingly complex problems are suitable for treatment by HarmonicBalance.jl, we consider in this section the harmonic equations for two coupled Duffing resonators (natural frequencies ω x and ω y , coupling k), governed bÿ In addition to obtaining steady states for multiple parameters using homotopy continuation, we solve the time dynamics to simulate an experimentally standard parametric sweep via an interface with DifferentialEquations.jl. This can be attained via the following lines of code @variables ω_x , ω_y , t , ω, F , γ , α , k , x ( t ) , y ( t ) ; free_eq = [ d ( d (x , t ) ,t ) + ω_x^2 * x + γ * d (x , t ) + α* x^3 -k *y , d ( d (y , t ) ,t ) + ω_y^2* y + γ * d (y , t ) + α* y^3 -k * x ] forces = [ F * cos (ω* t ) , F * cos (ω* t ) ] diff_eq = D i f f e r e n t i a l E quation ( free_eq -forces , [x , y ]) add_harmonic !( diff_eq , x , ω) # x will oscillate at ω add_harmonic !( diff_eq , y , ω) # y will oscillate at ω harmonic_eq = g e t _ h a r m o ni c_e qua ti ons ( diff_eq ) Execution of the above block also prints a detailed summary of the simplified system We visualise the amplitudes of stable steady-state solutions via plot ( soln , x = "ω" , y = " sqrt ( u1^2 + v1^2) " , plot_only =[ " physical " , " stable " ] , filename = " coupled_duffing_sols " ) ; Here we saved results to a .jld2 file by filling the filename keyword argument. Data can be subsequently recovered via use of the load function.
We are now in a position to check the accessibility of the steady states by performing a backward, time-dependent sweep of the frequency, initiating at a given particular condition: # select a solution and evolve from it s1 = soln [75] Steady-state solution amplitudes X 1 = u 2 1 + v 2 1 for a system initiated away from any steady state, rapidly collapsing into one of the available solutions. Adiabatic variation of the driving frequency follows a subset of stable, steady states. The arrows indicate the direction of the sweep; jumps occur at points where a solution branch ceases to be stable. The observed behaviour is hysteretic, i.e., it depends on the sweep direction. The parameters chosen are ω x = 1, ω y = 1.05, γ = 2 · 10 −3 , F = 10 −2 , α = 10 −1 , J = 5 · 10 −2 . The sweeping time is τ = 10 5 γ −1 , much slower than any relaxation time to minimise non-adiabatic effects.

Parametrically driven Duffing oscillator
Nonlinear effects in oscillating systems manifest not only when systems are driven through external forcing, as in Secs. 5.15.2 but also through periodic modulation of the system parameters, also known as parametric driving [77,101]. The simplest example of such an oscillator is the parametrically-driven Duffing oscillator or parametron [102][103][104][105][106], governed by the damped nonlinear Mathieu equation Interestingly, this system hosts finite oscillations when the parametric driving frequency approaches ω p ≈ ω 0 . In addition, above a critical driving rate, such a system is known to host two stable oscillation states with equal amplitudes but π-shifted phases. Their (subharmonic) response frequency, however, is at ω = ω p /2. We thus study the response at frequency ω via the following lines of code: @variables ω0 , γ ,λ,ωp ,Ψ ,α,η,t ,T , x ( t ) free_eq = d ( d (x , t ) ,t ) + γ * d (x , t ) + ω0^2*(1 -λ* cos (2*ω* t +Ψ) ) * x + α * x^3 + η * d (x , t ) * x^2 diff_eq = D i f f e r e n t i a l E q u ation ( free_eq , x ) add_harmonic !( diff_eq , x , ω) ; harmonic_eq = g e t _ h a r m o n i c_e qu ati ons ( diff_eq ) We can now study the 2D-phase diagram spanned by λ and ω (Fig. 6a). For it, we define the fixed parameters and specify the parameter range for λ and ω:  Note that such a representation reveals the transition from the region with a single stable solution to the region with two solutions corresponding to a pitchfork bifurcation, in the vicinity of which the imaginary part vanishes, while the real part splits [75].

Increasing computational complexity: Performance scaling
Here, we shortly illustrate the performance of HarmonicBalance.jl for varying system sizes. We consider a chain of linearly coupled Duffing oscillators, each with nonlinear damping (amplitude η, other parameters defined above Eq. (6)), x j (t) = F cos(ωt) , i = 1, 2, ..., N . (17) Similar systems have been explored in the context of combinatorial optimisation machines based on the mapping of effective spins to parametron networks [64][65][66][67][68][69]. The displacement of each oscillator, x i , is expanded using the harmonic oscillating at ω, giving the harmonic ansatz As explained in Section 2, each oscillator adds 2 harmonic equations of order 3. Therefore, a chain of length N and M = 1 leads to a set of 2N equations. The corresponding Bézout bound on the number of solutions, and hence the number of paths which must be tracked by the homotopy continuation algorithm, is 3 2N . In Table. 2, we show the computational times necessary for symbolic manipulation and homotopy solving as well as the number of complex and real solutions found. We highlight solving a chain of N = 5 resonators still takes a few minutes on a single CPU. Notably, the number of unique solutions is significantly smaller than the Bézout bound in all cases, with the number of real solutions being smaller still. The bottleneck is, therefore, the initial tracking of the homotopy from 3 2N paths to the relatively few non-singular ones. Subsequent steps such as solving for multiple parameter sets (i.e. tracking parameter homotopy paths) and determining stability (Jacobian evaluation and diagonalisation) only involve non-singular paths and thus are relatively inexpensive.
Two remarks regarding performance are in order. First, the process of path tracking is naturally well-suited for parallelisation, and HomotopyContinuation.jl does include the necessary functionality. Second, systems such as the Duffing chain possess spatial symmetries (in this case, inversion symmetry) and internal/gauge symmetries (e.g. discrete time translation, u i → −v i , v i → u i ) which cause some solutions to be degenerate. Making use of this property can further reduce the computational time needed.

Conclusion
In this work, we have introduced a Julia package to treat the steady-state behaviour of generic nonlinear, polynomial, dynamical systems with harmonic time dependence. Combining symbolic and numerical calculations with simple graphical capabilities, our package is developed with simplicity of use in mind while enabling the study of complex nonlinear systems.
The ability to explore the complete solution landscape with an ansatz that has a tunable level of complexity makes our tool ideal for working alongside experiments. Specifically, high-order processes (e.g., up-conversion) can be readily accounted for through the progressive addition of harmonic frequencies. We hope that the free availability of a user-friendly, high-performance code for harmonic nonlinear dynamics calculations will help advance multiple disciplines by allowing researchers to perform computations currently considered out of reach due to their complexity. In addition, we hope its implementation, sitting on Julia's rich ecosystem for high-performance scientific computing, will help nurture a multidisciplinary open-source community.

A Harmonic equations for a Duffing resonator
Here, we derive the harmonic equations for a single Duffing resonator, governed by the equation where we considered driving with a finite phase offset θ for generality. As explained in Section 2, for a periodic driving at frequency ω d and a weak nonlinearity α, we expect the response at frequency ω d to be dominant, followed by a response at 3ω d due to frequency conversion.

Expanding in ω d only
We first describe the established [4,77,82] approach to finding the steady states of Eq. (A.1), where frequency conversion is only added perturbatively. The starting point is a harmonic ansatz for x of the form Eq. (7), containing a single frequency ω d , with the harmonic variables u and v. The slow time T is, for now, equivalent to t. Substituting this ansatz into Eq. (A.1) results in ü + 2ω dv + u ω 2 0 − ω 2 d + 3α u 3 + uv 2 4 + F cos θ cos(ω d t) We see that x 3 in Eq. (A.1) generates terms that oscillate at 3ω d , describing the process of frequency upconversion. We now Fourier-transform both sides of Eq. (A.3) with respect to ω d to obtain the harmonic equations. This process is equivalent to extracting the respective coefficients of cos(ω d t) and sin(ω d t). Here the distinction between t and T becomes important: since the evolution of u(T ) and v(T ) is assumed to be slow, they are treated as constant for the purpose of Fourier transformation. Since we are interested in steady states, we drop the higher-order derivatives and rearrange the resulting equation to the form of Eq. (8) of the main text d d T Note that our assumption that u(T ) and v(T ) are slowly changing, i.e. composed of small frequency terms also sets constraints in Fourier space: we neglect all the frequencies which are not close to ω d . In the extreme case of constant u and v, the described frequency range reduces to the discrete frequency ω d . Steady states can now be found by setting the l.h.s. to zero, i.e., assuming u(T ) and v(T ) constant and neglecting any transient behaviour. This step is referred to in the literature as "balancing the harmonics" [79]. This results in a set of 2 nonlinear polynomial equations of order 3, for which the maximum number of solutions set by Bézout theorem is 3 2 = 9. Depending on the parameters, the number of real solutions is known to be between 1 and 3, see the solution diagrams Section 5. The steady states describe a response that may be recast as x 0 (t) = X 0 cos(ω d t + φ), where X 0 = u 2 + v 2 and φ = − atan(v/u). Frequency conversion from ω d to 3ω d can be found by setting x(t) ≡ x 0 (t) + δx(t) with |δx(t)| |x 0 (t)| and expanding Eq. (A.1) to first-order in δx(t). describes a simple harmonic oscillator, which is exactly soluble. Correspondingly, a response of δx(t) at frequency 3ω d is observed. Since this response is obtained perturbatively for each steady state of Eq. (A.4), no previously-unknown solutions are generated in the process.

Expanding in ω d and 3ω d
An approach in the spirit of harmonic balance is to use both harmonics ω d and 3ω d on the same footing, i.e., to insert the ansatz x(t) = u 1 (T ) cos(ω d t) + v 1 (T ) sin(ω d t) + u 2 (T ) cos(3ω d t) + v 2 (T ) sin(3ω d t) , (A. 6) with u 1 , u 2 , v 1 , v 2 being the harmonic variables. As before we substitute the ansatz into Eq. (A.1), drop second derivatives with respect to T and Fourier-transform both sides. Now, the respective coefficients correspond to cos(ω d t), sin(ω d t), cos(3ω d t) and sin(3ω d t). Rearranging, we obtain In contrast to Eqs. (A.4), we now have 4 equations of order 3, allowing up to 3 4 = 81 solutions (the number of unique real ones is again generally far smaller, as seen in Section 5). The larger number of solutions is explained by higher harmonics which cannot be captured perturbatively by the single-frequency ansatz. In particular, those where the 3ω d component is significant. Such solutions appear, e.g., for ω d ≈ ω 0 /3 where the generated 3ω d harmonic is close to the natural resonance frequency.

B Solving steady states for multiple parameter values
Evaluating Eq. (8) at each parameter set determines a distinct numerical system of polynomial equations, whose roots can be found using homotopy continuation; this generally requires tracking the maximum number of roots given by Bézout's theorem, which may be very timeconsuming. Having solved for one parameter set, the entire tracking procedure is not necessary for others -one only needs to track paths leading to non-singular solutions. In particular, we carry out 1. A "warm-up", i.e. path tracking from the roots of, e.g., U(u 0 ) = {u 2 −1, · · · } (d r , r = (1, 2, · · · , N M ) denotes the degree of the polynomialF r (u 0 )) towards the roots of a "generic" system, namelyF(u 0 ) evaluated with random complex (unphysical) parameters 10 . All singular paths are filtered out during this procedure.
2. The number of complex roots of the generic system is an upper bound for the number of roots ofF(u 0 ) = 0 for arbitrary parameters [30]. In a final step, we apply a parameter homotopy (a mapping between the collection of systemsF(u 0 ) = 0 obtained by only varying the parameters) towards a physical real parameter value. This may involve tracking a massively reduced number of paths, hence being much faster than step 1.
3. Repeat step 2 for each parameter value.
Despite the computational overhead of the problem's "complexification", complex solution paths allow for a stable path tracking algorithm, preventing singularities that will otherwise show up in a real-space parameter homotopy (see Ref. [30] for details).