The interior of neutron stars has been a subject of uncertainty for decades, as the interaction of supranuclear densities and relatively low temperatures is not currently replicable in terrestrial laboratories. In order to constrain the behavior of matter in such conditions, one typically turns to astrophysical observables of neutron stars, such as radius, mass, and tidal deformability. Most relevant to this thesis is the last quantity, which is inferred from the gravitational waves emitted by binary neutron star mergers. Current models of tidal interactions rely on quasi-static tidal deformations, which are only meaningfully approximate during the early part of the late inspiral. As the stars progress throughout the late inspiral stage, the tidal deformations become dynamical and nonlinear. This thesis aims to analytically model Newtonian tidal deformation contributions up to nonlinear terms. By identifying the derived dynamical corrections to the phase of the waveform, this thesis establishes a foundation for future modeling of linear and nonlinear PN contributions, ultimately aiming to help constrain the EoS of neutron stars using gravitational wave observations.

Since their first detection in 2015 by the LIGO, gravitational waves have revolutionized our understanding of compact objects. These dense stellar remnants such as neutron stars and black holes merge to generate gravitational waves. The merger of BNS can be accompanied by an electromagnetic signature of thermal emission called a "kilonova" at optical, near-infrared, and ultraviolet wavelengths. The BNS merger GW170817 led to the first detection of gravitational waves with electromagnetic radiation and the first confident detection of a kilonova. To date, only one multi-messenger BBH merger has been definitively detected. This contributes to large uncertainties in the relationship between binary parameters and properties of their kilonova counterparts, which are typically parameterized via fitting formulae. We attempt to place direct, data-driven constraints on the coefficients within the fitting formula through Bayesian analysis. We apply our framework to simulated BNS gravitational wave signals and their counterparts with a known model to place an independent constraint on the mapping between BNS parameters and kilonova properties. The constraint and this framework will demonstrate what will be possible with a population of many detected BNS mergers and counterparts in upcoming observing runs, given improved sensitivity.