## 1 Introduction  
The characteristic input-to-output behavior of analog audio circuits (timbre, transitory) rests on the possibly highly nonlinear components appearing in such systems. These components make the stability of the simulations difficult to guarantee. The motivation of this work stems from the following observations:  
1. Analog circuits combine energy-storing components, dissipative components, and sources. 2. Storage components do not produce energy, and dissipative components decrease it. In this sense, analog circuits can be considered as passive systems with external power supply. We shall exploit this passivity property by transposing it to the digital domain, ensuring the stability of the simulations (see [1; 2; 3]). The available approaches for the automated derivation of physical modeling and numerical simulation of audio circuits can be divided in two classes [4]: _wave scattering methods_ (WS) and _Kirchhoff's variables methods_ (KV). Mixed WS/KV methods have also been proposed in references [5; 6]. The well-established _wave-digital filter_ (WDF) formalism [7] belongs to the class of WS methods. For linear circuits, it provides a computationally realizable system of equations: First, by defining parametric wave variables for each elementary component and multiports (serial and parallel);Second, by discretizing the corresponding constitutive laws with the bilinear transform; and Third, by choosing the wave's parameters so as to reduce the computational complexity and to avoid instantaneous feedback loops. An extension to nonlinear circuits has been considered in [8; 9] and applied to, for example, the real-time simulation of vacuum-tube guitar amplifiers in [10; 11]. WDF ensures the passivity of the resulting digital system [7; 12], including systems with scalar nonlinearity [13]. However, the passivity property of WDF structures is not ensured for circuits with more than one nonlinear element (e.g., [14] SS3.2 and [15] S6). The class of KV methods for audio circuits encompasses nonlinear state-space representations [16]. Several modeling techniques are available to derive the discrete-time state-space model, either from the global time-continuous model (e.g., [17]) or from the interconnection of discretized elementary components (e.g., [18]). The resulting set of nonlinear implicit equations solved at each sample can be structured so as to obtain a computationally realizable system by applying the _K-method_ introduced in [19] with developments in [18]. However, this structure does not encode the passivity of the original circuit naturally, and it must be investigated on a case-by-case basis [20]. In this paper, we consider the _port-Hamiltonian systems_ (PHS) approach, introduced in the 1990's [21; 22; 23]. PHS are extensions of classical Hamiltonian systems [24], specifically defined to address open dynamical systems made of energy storage components, dissipative components, and some connection ports through which energy can transit. This approach leads to a state-space representation of physical systems structured according to energy flow, thus encoding the passivity property, even for nonlinear cases. This class of physical systems encompasses not only electrical circuits, but also multi-domain systems, such as loudspeakers, which involves electrical, magnetic, mechanical, acoustical, and thermodynamical phenomena. The port-Hamiltonian structure is derived by applying the Kirchhoff's laws to a given schematic, similarly to other existing approaches (e.g., WDF and K-method). Here, the advantage of the PHS formulation is the direct encoding of the underlying passive structure. This passivity property is transposed to the discrete-time domain by appropriate numerical methods, so as to ensure the numerical stability. For linear storage components (inductors and capacitors), the combination of the PHS structure with any of the trapezoidal rule or the mid-point rule yields the same numerical scheme that preserves the passivity in discrete-time.