## 1 Introduction
The frequency-based counterpart of Equation (1), constituting the so-called harmonic formulation, reads  
\[\left(-\omega^{2}\mathbf{M}+i\omega\mathbf{C}+\mathbf{K}\right)\mathcal{U}= \mathcal{F}, \tag{7}\]  
where \(\mathcal{U}\) and \(\mathcal{F}\) refer to the Fourier transform of the nodal displacement and forces, respectively, \(\mathbf{U}\) and \(\mathbf{F}\). As indicated, the main limitations of such a formulation are: (i) first, the treatment of nonlinear models that requires specific procedures; and second (ii) the necessity of solving the problem for each angular frequency \(\omega\) involved by the loading or present in the induced displacement (when nonlinearities are addressed). However, the use of harmonic analysis avoids the necessity of constructing the basis \(\mathbf{P}\) that needs the solution of the large-size eigenproblem related to  
\[\left(-\omega^{2}\mathbf{M}+\mathbf{K}\right)\boldsymbol{\phi}=\mathbf{0}, \tag{8}\]  
in which eigenvectors \(\boldsymbol{\phi}_{i}\) diagonalize the mass and stiffness matrix, as well as the one related to the damping as soon as a proportional damping is assumed. The main limitation of modal analysis is not only the necessity of solving such an eigenproblem for extracting the normal modes but the lack of validity of such basis in the case of parametric modelsin which the different matrices change when modifying the parameters grouped in vector \(\mathbf{\mu}\), i.e., \(\mathbf{M}(\mathbf{\mu})\) and \(\mathbf{K}(\mathbf{\mu})\). The solution of parametric eigenproblems remains a recurrent and still open problem. Another advantage of the harmonic formulation is the possibility of considering more general damping. In mechanics of materials, sometimes damping is assumed scaling with the inverse of frequency [4]. The interested reader can refer to Reference [5] that analyzed the theoretical consequences of assuming a frequency dependent dashpot parameter. In Reference [5], it was proved that, even if such a choice succeeded to fit experimental data, when coming back to the time-domain, causality is lost, and then the resulting expressions in the time-domain were called _non-equations_. Anyway, it is important to note that, even when considering complex nonlinear frequency dependent damping \(\mathbf{C}(\omega)\) in Equation (7), the problem in the frequency domain remains linear because, here, the frequency is a model parameter (or a model extra-coordinate within the Proper Generalized Decomposition (PGD) framework addressed later). In the context of model order reduction (MOR), in Reference [6], authors proposed a Proper Orthogonal Decomposition (POD)-based reduced order modeling operating in the time domain. [PERSON] and coworkers proposed an extension of their radial approximation [7] for addressing mid-frequency dynamics, the so called variational theory of complex rays [8]. In our former works, we considered a PGD (Proper Generalized Decomposition) formulation for constructing a parametric transfer function [9] that allowed efficient solutions of transient dynamics operating in the time-domain. On the other hand, the separation of variables, at the heart of PGD [10; 11], was extensively employed in the harmonic domain for solving multi-parametric dynamics, and was successfully extended to the non-linear case, and then combined with modal analysis [12]. In soil dynamics, the strain-stress relationship indicates that soil behavior during a seismic event is highly nonlinear and exhibits hysteresis [13; 14]. Nonlinear effects, are defined in terms of the shear modulus reduction and an increase of damping ratio as a function of the strain level. Various approaches are available to model this dependence, such as equivalent-linear (EL) or nonlinear (NL) analysis. Equivalent linear analysis is by far the most commonly utilized procedure in practice, however, may lead to a conservative approach [15].