Tidal Wave Breaking in the Eccentric Lead-In to Mass Transfer and Common Envelope Phases

The evolution of many close binary and multiple star systems is defined by phases of mass exchange and interaction. As these systems evolve into contact, tidal dissipation is not always sufficient to bring them into circular, synchronous orbits. In these cases, encounters of increasing strength occur while the orbit remains eccentric. This paper focuses on the outcomes of close tidal passages in eccentric orbits. Close eccentric passages excite dynamical oscillations about the stars' equilibrium configurations. These tidal oscillations arise from the transfer of orbital energy into oscillation mode energy. When these oscillations reach sufficient amplitude, they break near the stellar surface. The surface wave-breaking layer forms a shock-heated atmosphere that surrounds the object. The continuing oscillations in the star's interior launch shocks that dissipate into the this atmosphere, damping the tidal oscillations. We show that the rapid, nonlinear dissipation associated with the wave breaking of fundamental oscillation modes therefore comes with coupled mass loss to the wave breaking atmosphere. The mass ratio is an important characteristic that defines the relative importance of mass loss and energy dissipation and therefore determines the fate of systems evolving under the influence of nonlinear dissipation. The outcome can be rapid tidal circularization ($q\ll1$) or runaway mass exchange ($q\gg1$).


INTRODUCTION
In binary and multiple star systems, the combined evolution of the stars and the orbits they share can drive systems that begin separate to directly interact. The resulting interactions include stable and unstable mass transfer, common envelope phases in which one object engulfs its companion, and complete mergers of two or more objects. Indeed, these sorts of interactions define the evolution of more than half of massive stars (Sana et al. 2012;Duchêne & Kraus 2013;Moe & Di Stefano 2017). In order to model how populations of multiple stars evolve and produce remnants, we need predictive theories for the mapping systems through these phases. With so many possible outcomes, performing and validating these theoretical models has been challenging and many open questions remain (e.g. Paczyński 1971;Hur-ley et al. 2002;Ivanova et al. 2013;Pavlovskii & Ivanova 2015;Broekgaarden et al. 2021).
The properties that describe the initial conditions for the mapping of multiple-star systems through phases of interaction are 1) the initial state of the system orbit, and 2) the structural properties of the objects themselves. Together, these properties define the outcome of a subsequent hydrodynamic interaction. In this paper, we focus on the property of initial orbital eccentricity at the outset of mass exchange. Pairs of stars within multiple systems are formed with a wide distribution of initial orbital eccentricities. These initial distributions are modified by processes of tidal dissipation into the stellar interiors, gravitational torques from other bodies (e.g. secular torques in hierarchical multiples), and impulsive mass loss phases like supernovae.
For any of these reasons, a system might evolve toward a phase of mass exchange while in an eccentric orbit (e.g. Dosopoulou & Kalogera 2016a,b). In these situations, tidal dissipation competes with that evolutionary process to determine whether a mass exchange interac-tion begins in a circular or eccentric orbit (Verbunt & Phinney 1995;Hurley et al. 2002). Vigna-Gómez et al. (2020) identified, for example, that in pre-commonenvelope phases involving systems that go on to produce merging double neutron star systems, the characteristic timescales for tidal circularization and spin synchronization are often longer than the radius growth timescale of the giant star in the system, apparently driving the system toward a common envelope phase while still in an eccentric orbit (Staff et al. 2016;Shishkin & Soker 2020;Glanz & Perets 2021;Schreier et al. 2021). Vick et al. (2021) analyzed this scenario with more sophisticated models of tidal dissipation into convective stellar envelopes (Vick & Lai 2020) and found that the action of tides circularizes some systems, but many remain eccentric and asynchronous by the time the giant star evolves to overflow its Roche lobe at periapse. Another representative example of this class of interaction comes when secular torques drive a system toward Roche lobe overflow in an eccentric orbit, in this case tidal dissipation competes with the torque on the orbit to determine the orbital properties (Naoz & Fabrycky 2014;Naoz 2016;Salas et al. 2019).
As a system evolves toward the extreme of mass transfer in an eccentric or asynchronous orbit, the primary star experiences dynamical tides of increasing strength. These tides are dynamical because they cannot be in equilibrium with the time-varying tidal potential as viewed from the stellar fluid frame. In this paper, we show that as these tides increase in amplitude and waves associated with the tidal oscillations break. In even closer passages, the amplitude becomes so large that mass is transferred from the primary to the perturber. Wave breaking dissipates the organized motion of the oscillation into random motions, through a turbulent cascade, heat.
Wave breaking, therefore, dissipates the energy deposited into dynamical tides when the oscillations reach an amplitude that leads them to become nonlinear. We perform hydrodynamic models of the excitation of fundamental oscillation modes by an eccentric tidal passage, and demonstrate the wave breaking that occurs when these waves reach unsustainable amplitudes. We extend our models for many dynamical times of the primary star after the periapse passage in order to quantify the rate and spatial distribution of energy dissipation associated with the damping of the fundamental mode by wave breaking. Through this dissipation, wave breaking shapes the crucial initial conditions of system orbits and the structures of tidally-excited stars at the onset of Roche lobe overflow for eccentric systems approaching mass exchange.
In Section 2, we describe our hydrodynamic simulation methodology and the linear theory that forms a baseline for interpreting our numerical results. In Section 3, we describe the process of tidal wave breaking and the dissipation of mode energy into turbulence. In Section 4, we discuss some potential implications of our results, and in Section 5, we conclude.

Hydrodynamic Models
The simulations presented in this work are performed using the Athena++ software ). Our method is based on that of MacLeod et al. (2018aMacLeod et al. ( ,b, 2019; MacLeod & Loeb (2020a,b). We model the interaction of a tidally-perturbed stellar envelope with two point masses, one representing the excised center of the star, and the other representing the perturbing object. Our simulations are three dimensional, and are conducted in spherical polar coordinates with their origin at the center of the star.
The Eulerian approach of our Athena++ simulations allows high resolution in the low-mass regions at the limb of the star. This high resolution is important because it allows us to capture the shock-heating that results from wave breaking and the formation of the lowmass, high entropy atmosphere around the star. Though there are many advantages to Lagrangian hydrodynamics in stellar modeling, these mass-based methods do not have the ability to resolve the crucial features that we study in this work.
The total mass of the star is M 1 = m 1 +m g , where m 1 is the excised core mass, and m g is the gaseous envelope mass. The perturber has mass M 2 , and the mass ratio is q = M 2 /M 1 = 0.1. The original radius of the star is R 1 . Our calculations are performed in units where G = M 1 = R 1 = 1, but may be rescaled to any characteristic stellar M 1 and R 1 . The unit velocity is v 1 = GM 1 /R 1 and the unit time is t 1 = R 3 1 /GM 1 . The initial condition of the stellar envelope is a nonrotating spherically-symmetric polytropic envelope surrounding the excised core, which has mass m 1 = 0.1M 1 . This envelope is constructed by choosing a pressure and density at 0.1R 1 , then integrating outward to R 1 following the equations of of hydrostatic equilibrium closed by the polytropic index P ∝ ρ Γs , where Γ s = γ = 1.35 or 5/3, thus implying that the envelope is isentropic in either case. The initial density and pressure are iteratively chosen such that M 1 = 1 at R 1 = 1. Figure 1 shows the structure of our model objects. Qualitatively, these models represent giant-like stars with central condensations of mass and isentropic envelopes with differing degrees of central condensation. Outside of R 1 , the solu- . Structure of two model polytropic stars used in our simulations. The upper panel shows density, the central enclosed mass, and the lower the radial wavenumber kr ≈ ωα/cs, adopting the frequency of the l = 2 fundamental mode. Both objects have a core mass of 0.1M1 enclosed within 0.1R1, surrounded by polytropic envelopes with structural indices of Γs = 1.35 and 5/3. The Γs = 1.35 case is more centrally-concentrated, with less mass in its outer envelope while the Γs = 5/3 case is more homogeneous.
tion is joined with a low-density hydrostatic background of constant sound speed of v 1 /3.
Our frame of reference orbits with m 1 and is therefore accelerated, but is non-rotating. The coupled equations of inviscid gas dynamics that we solve are expressing mass continuity, the evolution of gas momenta, and the evolution of gas energies. In the above expressions, ρ is the mass density, ρv is the momentum density, and E = + ρv · v/2 is the total energy density and is the internal energy density. Additionally, P is the pressure, I is the identity tensor, and a ext is an external acceleration (source term). The equations above are closed by an ideal gas equation of state, P = (γ − 1) , in which γ is the gas adiabatic index or ratio of specific heats. The acceleration source term on the right-hand side of equations (1), a ext , represents the forces associated with the binary, as well as fictitious accelerations associated with our choice to perform the calculation in a noninertial frame of reference. The source term contains the gravitational accelerations from point masses m 1 and m 2 . The acceleration of m 2 is softened across 0.1R 1 , using a spline kernel (Hernquist & Katz 1989). It also contains the self-gravitational acceleration of the gas on itself. This we model via a monopole approximation by summing the enclosed mass within a given radius from m 1 , and setting the self-gravitational term acceleration equal to −Gm g (< r)/r 2 . Finally, we include the inverse of the acceleration of m 1 to compensate for our choice of non-inertial reference frame. A more detailed account of these terms is given in equations (7) Our mesh extends over the full 4π steradians of solid angle and from 0.1R 1 to 100R 1 . The zones are uniformly spaced in θ and φ coordinates but logarithmically spaced in r so that zone shapes remain roughly cubic across a wide range of radial scales. The base mesh is made up of 480x192x384 zones in r, θ, and φ. A nested layer of static mesh refinement increases the resolution by a factor of two between 0.7R 1 and 1.3R 1 . This mesh is divided into blocks of 32 3 zones for parallelization.
The inner radial boundary surrounding m 1 is located at 0.1R 1 , and is reflecting. The outer radial boundary at 100R 1 allows outflow from the mesh but not inflow onto the mesh. The boundaries at θ = 0 and θ = π are polar boundaries that allow material to flow through the coordinate singularity at the poles. The boundary at φ = 0 and φ = 2π is periodic, to represent the continuous flow across this coordinate transition. Near the θ = 0 and θ = π coordinate singularities we average conserved quantities across multiple zones in the φ direction to create effectively larger, more cubical zones that avoid extreme aspect ratios near the poles. This averaging preserves the data structure of the mesh while preventing excessive φ resolution, which can severely restrict the courant-limit on the time step, in the immediate vicinity of the polar boundary.
The orbit of the perturber and star are initialized such that the perturber follows an eccentric orbit with a = 100 and r p = a(1 − e) in the range of 1.1R 1 to 1.9R 1 . We integrate the orbit from apoapse to a separation of 10R 1 treating both masses as point masses, then initialize the hydrodynamic model. Particle-gas interactions are computed by direct summation over the simulation zones, and time integration is done by a leapfrog method with timestep equal to the hydrodynamic timestep. As the model is initialized, we relax the hydrostatic profile onto the mesh, progressively turning off a damping term for 5t 1 , then turn on m 2 over the subsequent 0.1t 1 . The orbital integration then begins and the star proceeds toward its periapse passage with the perturber. At our fiducial resolution, these choices allow us sufficient fidelity to capture exchanges of energy and angular momentum with dimensionless amplitudes 10 −6 . Figure 2 shows the exchange of angular momentum (top panel) and energy (bottom panel) between the stellar fluid and the orbit in our model. We plot this exchange in the weakest encounter for Γ s = 1.35, which has the smallest amplitude exchange of the models we report on. Because conservation is not guaranteed in our algorithm, the departure of these quantities from perfect conservation traces the accumulation of error in our algorithm. The overall conservation of system angular momentum in our simulations is less well preserved than the exchange, with accumulated error in the overall integra-tion of the orbital motion accumulating into an angular momentum error at the ∼ 10 −4 level after 100t 1 .
We incorporate scalar fields that are advected with the flow in order to provide Lagrangian diagnostics of the changes in energy and momenta imparted to particular fluid parcels. At the end of the relaxation period of damping the initial profile onto the mesh, we assign scalar tracers to represent material originating within the stellar envelope (to distinguish from the numerical background material). We also trace the initial r 0 , θ 0 , φ 0 , coordinates, as well as the initial potential, kinetic and internal energies. These tracers then passively follow the flow, allowing us to measure Lagrangian displacements of fluid in either physical or energy space despite our Eulerian method of solving the fluid equations.

Linear Theory
We compare our numerical results with predictions for the radial velocity field from linear tidal theory (see Figure 4). In the linear theory, the response of the object to a tidal potential is a sum of the response of each eigenmode. Eigenmodes are indexed α = {n r , l, m}. We use the stellar oscillation code GYRE to calculate the l = 2 − 10 f-modes (n r = 0) of the stellar model (Townsend & Teitler 2013). The tidal energy transfer in a pericenter passage, ∆E α , is given by equation (21) of Vick et al. (2019) assuming a non-rotating object. This calculation is similar to those of Press & Teukolsky (1977) and Lee & Ostriker (1986) but allows for an eccentric orbit rather than a parabolic encounter. We compute the mode masses M α using equation (9) of MacLeod et al. (2019). Equation (3) of that paper relates ∆E α and M α to the radial velocity perturbation δv α (R 1 ). The predicted power spectrum from the linear theory is then given by the sum where the factor of (4π) −1 arises from the normalization of the spherical harmonics. Important points of reference for comparison of Lagrangian perturbations the oscillating stellar fluid include the displacement, ξ, relative to the local mode wavenumber, k, and the velocity, δv, relative to the local sound speed. The dispersion relation of plane acoustic waves is k 2 = ω 2 /c 2 s , which relates k to the oscillation frequency and local sound speed. We will focus near the stellar surface, where the radial component of k is most important. We make the order of magnitude approximation that we can treat the mode near the surface as a plane wave, for which Appendix E of Aerts et al. (2010) includes a more complete derivation of the local wavenumber in the asymptotic regime where the wave properties vary rapidly against the background, which is satisfied near the stellar surface. We find that equation (3) is an accurate approximation in the region of interest near the stellar limb, and adopt this simplification in what follows. The value of k r , approximated by equation (3), is plotted for our model objects in Figure 1. Although the fundamental mode has no nodes, near the stellar surface k r R −1 1 .

WAVE BREAKING IN CLOSE, ECCENTRIC PASSAGES
In a close tidal passage, the time-dependent gravitational potential of the perturber excites tidal oscillations in the stellar fluid. In the extreme of this case, fundamental modes excited by the tide break -implying the irreversible deformation of the wavelike motion of the oscillation -as the perturbing force is removed and the tidal wave peak crashes back to the stellar surface. As we discuss, the result of this irreversible deformation is dissipation. Figure 3 shows one such periapse passage with tidal wave breaking. In this case, a perturber of M 2 = 0.1M 1 passes at a minimum distance r p = 1.3R 1 . For context, this system would overflow its Roche lobe in a synchronous, circular orbit at a ≈ 1.7R 1 (In an eccentric, asynchronous orbit the critical separation for Roche lobe overflow is smaller Sepinsky et al. 2007). In the isentropic envelope, modes without radial nodes, the fundamental modes of various degree and azimuthal order, l, m, are the primary modes excited by a tidal passage. These modes have eigenfunctions with peak amplitude at the object's surface.

Periapse Passage
The amplitude of the tidal wave that is excited during the periapse passage is so large that the associated radial velocities of the oscillation are greater than the local gas sound speed. With this supersonic motion, the wave loses phase coherence causing material to self-intersect and crash back to the star's surface as the perturbing force is reduced. This fallback launches shocks through the star's outer layers, highlighted by purple regions of strong velocity convergence in the center panels of Figure 3. Their effect is to increase the specific entropy of the shocked material, traced by ln(P/ρ γ ) in the lower panels.
We examine the process of fundamental mode wave breaking in this same r p = 1.3R 1 model in more detail by analyzing the star's fluid motion in a spherical surface at it's original radius r = R 1 . Figure 4 decomposes the radial velocity field at r = R 1 into its spherical harmonic components. We plot the power in each angular order l summed over azimuthal orders m, where Figure 4 shows a transformation of the power spectrum near the time of periapse, t p . The initial mode spectrum, as observed from (t − t p )/t 1 = −2 to 0, mirrors the spectral shape predicted by the linear theory, with increasing amplitude as periapse approaches. At (t − t p )/t 1 > 1, the spectral shape has transformed. The peak still lies in the l = 2 quadrupole mode, but the power in higher l falls off roughly as a power law. We find that the l −5/3 scaling predicted by Kolmogorov turbulence (Kolmogorov 1941) provides an approximate description of the mode spectrum for l > 2. We caution, however, that our preliminary investigations have found that this power law is not universal. For example, in cases of much larger mass ratio, where the tide is very symmetric, we see a pronounced even-odd pattern in the mode power spectrum, with substantially more power in the even modes, even after wave breaking. This transformation of the mode power spectrum by wave breaking represents the extraction of coherent energy in lower degree modes (l = 2 and l = 3) and its deposition into disordered motion that adds power to higher l modes. Before periapse, the energies in all of the modes are increasing, with their relative powers remaining nearly constant. Near the periapse passage, while the l = 2 and l = 3 power continue to rise, the l > 4 powers momentarily flatten then dramatically increase near (t − t p )/t 1 = 2, as the tidal wave begins to break and shocks are launched through the material, modifying the wavelike radial velocity pattern. The slices at r = R 1 of radial velocity, labeled a, b, and c in Figure  4, visualize the smoothly varying wave pattern in panel (a), the beginning of wave breaking and shocking in the abrupt transition in panel (b), and the transfer of power smaller-scale velocity features, or higher l, m modes, in panel (c).

Conditions for Wave Breaking
. Wave breaking during a close periapse passage. These panels show slices through the equatorial plane. Top panels show gas density, center panels show velocity divergence, and lower panels show specific entropy. Contours in each panel are at ρ = 10 −5 , 10 −4 , and 10 −3 M1/R 3 1 . The dashed line and cross marker indicate the orbital path and position of the perturbing object, respectively. In this case, the perturber passes at a periapse distance of 1.3R1, giving rise to a large amplitude tidal wave, which collapses on itself and breaks, forming shocks. These shocks are traced by purple regions of strong velocity convergence. Shocks increase the entropy of some of the initially isentropic stellar material. Figure 5 explores the conditions for tidal wave breaking of the fundamental l = 2 mode in simulations with different periapse distances from r p = 1.1R 1 to r p = 1.9R 1 . The left panels show the distortion of the star's surface by the passage of the tidal perturber. Near periapse (left panel) we see a tidally-excited wave, lagging the position of the perturbing body, that increases in amplitude the closer the periapse distance. In the second panel, taken when the orbital phase has advanced slightly, we see that the surface contours in the r p = 1.1, 1.3, 1.5, and 1.7R 1 cases are no longer smoothly wave like, but instead host sharp reversals and multiple peaks. By comparison to the fluid properties seen in Figure 3, we see that these surface discontinuities are the signatures of wave breaking. By contrast, the contours of the r p = 1.9R 1 case remain smoothly wavelike. For this particular scenario, we therefore see wave breaking occurring at periapse distances less than the separation where a synchronously-rotating star would overflow its Roche lobe in a circular orbit.
In the right-hand panels of Figure 5, we look at the mass distributions of stellar material in radial Lagrangian perturbation, ξ r , measured in snapshots that are within a half dynamical time of periapse, 0 < (t − t p )/t 1 < 0.5. Closer periapse passages yield larger stellar distortions and ξ r . As described in Section 2.2, a point of comparison for ξ r is the radial wave number, k r ∼ ω α /c s , where the approximation is valid near the star's surface (Aerts et al. 2010). The condition |kξ| > 1 is often written as an approximate criteria for wave breaking, implying that waves tend to break when their amplitudes become large relative to their scale length (Barker & Ogilvie 2010;Barker 2011;Su et al. 2020). In our case, wave breaking is driven by the large radial displacements of the tide near the surface, so the condition becomes |ξ r | 1/k r or |k r ξ r | 1.    . The transformation of the spherical harmonic mode spectrum during the periapse passage due to wave breaking, in the rp = 1.3R1 case. The upper panel compares the mode spectrum to that predicted by linear theory. In the lead-in to periapse, the simulated mode spectrum mirrors the linear prediction. In the immediate aftermath the mode spectrum is transformed to a power law, similar to l −5/3 , by wave breaking. In the lower panel, we see how, as the power in low-degree modes grows, wave breaking transforms the radial velocity field at r = R1, extracting power from the coherent, wave-like motion and depositing it in higher-degree modes and turbulence.
The mode velocities implied are |δv r | = ω α |ξ r |, which means that the the wave breaking condition can be reexpressed as a condition on supersonic radial mode velocities, |k r ξ r | ≈ |δv r |/c s .
The right hand panels of Figure 5 compare ξ r to the inverse surface radial wave number, 1/k r (R 1 ), and examine the distribution of |k r ξ r |. We see that in each of the cases where we note the presence of wave breaking (r p ≤ 1.7R 1 ), there is material with ξ r > 1/k r (R 1 ). Further, the distributions of |k r ξ r | show that the amount of stellar mass exceeding the breaking condition of |k r ξ r | > 1 increases in stronger periapse passages. By contrast, in the r p = 1.9R 1 case without wave breaking, |k r ξ r | < 1 for all of the star's material. It is, however, worth noting that even when some material exceeds the breaking condition, the vast majority of the stellar material remains in the linear regime where |k r ξ r | 1. For example, in the strongest encounter we simulate (r p = 1.1R 1 ), only ∼ 10 −3 M 1 is involved in the wave breaking.

Shock-Heated Atmosphere
Just as wave breaking leads to dissipation in ocean waves, wave breaking on the surface of the star limits the maximum sustainable amplitude, |ξ r | 1/k r (R 1 ) (or, equivelently, |δv r | c s (R 1 )), and dissipates tidal oscillations. The wave breaking layer damps the still-coherent linear oscillations deeper in the star's interior. The snapshots of Figure 6 visualize the r p = 1.3R 1 model across the 80 dynamical times following its periapse passage. In the leftmost panels, at 10 dynamical times after periapse passage, there are still abundant waves on the object's surface, and the surrounding atmosphere is quite chaotic. From the velocity divergence, we see that oscillations in the object's interior drive shock waves into the atmosphere. These form spiral patterns as they propagate both radially outward and azimuthally around the object with the pattern speed of the mode. As these shocks dissipate into the atmosphere, it develops a pattern of coherent rotation at v φ ∼ 0.4v 1 .
The spherical averages of Figure 7 exhibit the development of the rotating atmosphere layer surrounding the star, from snapshots approximately 80 dynamical times after the periapse passage. Dotted lines show the profile of the initial conditions of our simulations, which include the physical object interior and an artificial hydrostatic atmosphere that is included because vacuum conditions are not permitted in our hydrodynamic method. In each of the periapse passages with r p ≤ 1.7R 1 , an atmosphere that is significantly denser than the initial condition develops overlying a relatively unmodified interior of the star. Because this layer is shock-heated, it has specific entropies higher than the initially-isentropic interior of the star. The atmospheres also have similar profiles of Mach number ( M 1) and azimuthal velocity (v φ ∼ 0.4v 1 ). Meanwhile, the interior of the star maintains nearly its initial profile of density, constant specific entropy, and no coherent rotation.
Shock waves launched by the nodes of the oscillatory waves in the star, as visualized by Figure 6, create these shared properties of atmospheric mach number and rotation. The l = 2 fundamental mode propagates around the star at an azimuthal velocity v φ = ω 2 R 1 ≈ 2.26v 1 . This represents the azimuthal velocity of the shock. Taking the strong shock limit, the post shock velocity would be at least v φ = ω 2 R 1 (γ − 1)/(γ + 1) ≈ 0.34v 1 . ] r p = 1.9R 1 r p = 1.7R 1 r p = 1.5R 1 r p = 1.3R 1 r p = 1.1R 1 Figure 5. Contours indicating the stellar surface (left panels) and mass distributions with respect to Lagrangian displacement ξr (right panels). Through more distorted surface contours and larger ξr, we observe that closer tidal passages lead to larger amplitude tidal waves near the star's surface. In the rp = 1.1R1 through rp = 1.7 cases, this wave breaks, and we see abrupt distortions of the wavelike profile of the surface contours. The weakest, rp = 1.9R1 case features a sufficiently small wave amplitude that it does not break, and the object surface retains its oscillatory motion. From the mass distributions, we see that the condition |krξr| > 1 distinguishes the wave breaking cases from those that remain linear.
Thus, the development of rotation in the atmosphere is related to the shock heating process and is a direct consequence of the dissipative action of shocks launched by the oscillation. Although we have focused on the spherically-averaged properties for simplicity, this net rotation represents angular momentum in the atmosphere, which forms into a thick torus (scale height divided by radius of h/r ∼ 1). The mass of the atmosphere layer can be measured by several criteria involving the transformation of properties from the interior to the atmosphere. In the righthand panel of Figure 7, we compare the masses exceeding the wave breaking criterion to the mass that forms the emergent atmosphere. At late times, approximately 80 dynamical times after periapse passage, we measure the shock-heated and supersonic (M > 1) atmosphere masses. The shock heated masses are determined by a specific entropy jump equivalent to a shock of M ≥ 2. These atmosphere masses are compared to masses of material that exceed wave breaking criteria at periapse, as visualized by the mass distributions of Figure 5. We find that the non-linear mass, |k r ξ r | > 1, at periapse is very similar to the supersonic atmosphere mass at late times. The linear theory prediction of the mass with |k r ξ r | > 1 is similar, but lower by a factor of a few, and does not predict wave breaking in the r p = 1.7R 1 case. This difference represents the importance of nonlinearity in these strong tidal passages.
Meanwhile, the total shock heated atmosphere mass at late times is more closely traced by |k r (R 1 )ξ r | > 1 at periapse, which adopts the wave number at the star's surface. This criterion implies |δv r | c s (R 1 ), or that oscillatory motion is more rapid than the surface sound speed. The qualitative explanation appears to be that when tidal waves break and fall back to the surface, they shock the layers of material with sound speed lower than their characteristic velocity.  Figure 6. Damping of oscillatory motions by dissipation into the turbulent atmosphere. While wave breaking during the periapse passage initially increases the specific entropy of the surface layer, the subsequent interaction of the primary l = 2 oscillatory mode in the interior with this overlying atmosphere drive the damping of the oscillation, and the development of a coherently rotating surface layer subject to spiral shock waves that propagate azimuthally with the l = 2 mode pattern speed. Figure 8 shows the power in the l = 2 modes over the simulation durations, evaluated at the original donor radius, r = R 1 . After reaching a peak shortly after periapse, the mode power displays some variations about an overall decrease over time. When we measure the mode power at other locations in the star's interior, we see that it changes in proportion. Thus, the decrease over time represents a dissipation of the mode's energy, rather than a local change at r = R 1 . The variations come from exchanges between modes. One such example can be seen clearly in Figure 4, where at t − t p ≈ 6t 1 the l = 3 power increases as the l = 2 power decreases, then they reverse.

Dissipation
The deepest periapse passages, which experience the most extreme wave breaking, undergo the most rapid dissipation, even once their oscillatory amplitude has decreased. This means that the heirarchy of mode power by the end of the simulations is reversed. The initially strongest tides of the r p = 1.1R 1 model have the smallest residual amplitude after 80 dynamical times. The dissipation is exponential, following a constant decay time, τ diss . In Figure 8, we overplot best fit exponential models for the mode power for t − t p > 2t 1 . In these fits, S f f (l = 2) ∝ exp [−(t − t p )γ nl ]. The right hand panel of Figure 8 shows that γ nl is largest (most rapid dissipation) for r p = 1.1R 1 and becomes progressively smaller for the larger periapse distances, following ] periapse, |k r r | > 1 periapse, |k r (R 1 ) r | > 1 late time, shock heated late time, > 1 linear theory, |k r r | > 1 Figure 7. Spherically-averaged profiles at approximately 80 dynamical times after periapse passage (left panels). We see that wave breaking has disturbed the outer layer of the object's profile adding kinetic energy (which has both rotational and turbulent components) and heat, as traced by the specific entropy. The right panel shows the mass of the late-time atmosphere compared to the non-linear mass by several measures. We find that the linear theory prediction and simulated periapse mass with |krξr| > 1 correlate with the late time mass with M > 1. The total late time shock-heated mass correlates with the periapse mass with |kr(R1)ξr| > 1, which implies |δvr| > cs(R1).  simulation periapse, |k r r | > 1 linear theory, |k r r | > 1 Figure 8. Power in the l = 2 mode as a function of time and its damping. After periapse passage, we see that the mode decays exponentially, with characteristic timescale that is shortest in the strongest encounters. This effect is so extreme that after ∼ 40 dynamical times, the closest periapse passage has weaker ongoing oscillations than the furthest periapse passage. These strong encounters generate shorter dissipation timescales due to the larger masses of shocked atmosphere material generated by the initial wave breaking. The second panel shows the dissipation rate associated with the exponential fits (the rp = 1.9R1 case is shown with a lower limit because of background material in our hydrodynamic models). The overall behavior is similar to that of our estimate of the dissipation efficiency, with th ∼ 1 and matm equal to the wave-breaking mass at periapse (|krξr| > 1).
a roughly exponential trend with periapse distance. We note here that our measurement for r p = 1.9R 1 is numerically dominated, in that it comes primarily from the background material in the initial condition of our hydrodynamic models rather than material generated by wave breaking (see the angle-averaged density in Figure  7). Understanding how shocks dissipate energy from the oscillation into the atmosphere can lead us to a prediction for how the oscillatory modes dissipate energy over time. The oscillation energy of a fundamental mode indexed by α = {l, m} is defined as where M α is the mode mass, an integral quantity of the mode eigenfunction and stellar density profile, and δv r (R 1 ) is the characteristic velocity amplitude of the mode at the star's equator. Because the mode amplitude and also the torus density decrease toward the poles, dissipation is equatorially concentrated. However, we might expect a prefactor of the order of unity to represent this geometric arrangement. Therefore, in cases where wave breaking occurs, if shock waves dissipate a fraction th ∼ 1 of the specific energy of δv r (R 1 ) 2 /2, we can estimateĖ whereṁ atm is the rate that shocks sweep through atmosphere mass, m is the azimuthal order of the mode (and thus the number of shocks launched around the circumference of the star), and Ω atm is the atmosphere rotation frequency. The characteristic decay rate of the mode due to this nonlinear dissipation, γ nl , can then be estimated as In our case, the most excited modes are the l = 2, m = ±2 fundamental modes, and Ω atm 0.3, so we estimate γ nl ∼ 0.15 th (m atm /M α ). We use our estimates of the wave-breaking mass at periapse, with |k r ξ r | > 1, as reported in Figure 5, and M α (l = 2) ≈ 5.5 × 10 −4 M 1 , to estimate dissipation timescales. We find that a factor th ≈ 1 provides a good fit for the normalization and overall trend of our models with differing periapse distances. Thus, damping of mode energy comes about as waves that propagate freely in the stellar interior convert to dissipative  shock waves at the atmosphere boundary. The relative masses of the mode, M α , and the dissipative atmosphere, m atm ∼ m(|k r ξ r | peri > 1), determine the damping rate. The predicted atmosphere masses from linear theory also exhibit a similar trend and normalization to the simulation models, but do not capture the presence of wave breaking in the r p = 1.7R 1 case.

Dependence on Stellar Structure
In this section, we briefly examine the sensitivity of our conclusions to the structure of the star, using its polytropic index as a proxy. We compare the Γ s = 1.35 models that have been the basis for our discussion so far to an otherwise-identical polytrope model with Γ s = γ = 5/3. This less-compressible polytrope has a less centrally concentrated density profile, and serves as a comparison to study how object structure affects the energy deposition and wave breaking mass.
We use our simulated periapse passages to measure the energy exchange between the mode and the tide along with the wave-breaking mass that results as a function of a given mode energy. These quantities are shown in Figure 9. In the upper panel, we show the energy deposited in a single periapse passage, ∆E, and wave-breaking mass m wb , as a function of periapse distance for each of our simulated objects. In this panel m wb is the mass that has |k r ξ r | > 1 at periapse (See Figure 7). We compare (in dotted lines of matching color) to the linear theory predictions, summing modes from l = 2 to l = 10. At a given periapse distance, more energy is deposited into oscillations of the Γ s = 5/3 models, because these less centrally-concentrated structures have more mass in the outer layers that are most perturbed by the tidal potential. In both cases, the linear theory provides an accurate description of ∆E(r p ). By contrast, the linear theory underpredicts the mass in the wave-breaking layer, because this mass arises from the high-energy tail of non-linearly perturbed material, as shown in the periapse distributions of Figure 5. This is particularly evident in the more mass-rich outer layers of the Γ s = 5/3 model.
The lower panel of Figure 9 shows the mass in the wave breaking layer as a function of energy deposited into tides. Interestingly, both polytropic structures arrive at similar results, in which we observe a roughly linear dependence between the wave breaking mass and mode energy, above a certain non-linearity threshold. This linear dependence can be cast in terms of a specific energy, as shown by the background contours in Figure 9, which shows that the ratio of energy deposited over wave breaking mass is between the GM 1 /R 1 and 10GM 1 /R 1 . This, however, does not imply that the specific energy of material in the wave breaking layer is that high, in fact it is usually at least one to two orders of magnitude lower because most of ∆E is deposited in the more mass rich lower layers at specific energies in the linear regime, and the mass exceeding the wave breaking condition represents the tail of the energy distribution.
These results imply that there may be a relationship between energy deposited into the tidal oscillation and mass in the wave-breaking layer that does not depend on stellar structure. In the lower panel of Figure 9 we contextualizes this finding with a linear approximating form, where C wb = 0.25 is a constant, and E nl ≈ 2 × 10 −5 GM 2 1 /R 1 is the mode saturation energy that im-plies wave breaking at higher energies. The upper limits in Figure 9 represent estimated mass-resolution limits for our simulations in which we do not observe wave breaking. In principle unresolved wave breaking could occur if the outer scale heights of the object were more finely resolved. In practice, there is likely some transition around the branches of the approximation above. An exponential cutoff of the wavebreaking mass is a physically-likely case due to the exponential density profile of the outer layers of the stellar models (and the behavoir of the linear theory predictions). However, it is not possible to constrain the shape of this cutoff with our current models, so we focus on the two clearly defined regimes of equation (10).

Comparison to Linear Tidal Dissipation: Location and Magnitude
How tidal dissipation occurs, and precisely where the dissipation is located in a tidally-excited object has long been seen as a crucial question in the fate of tidallycircularizing objects (Fabian et al. 1975;Podsiadlowski 1996). For many objects, such as giant planets around their host stars, or stars around more-massive black hole companions, the tidal object must mediate more than its own binding energy in orbital and oscillation energy exchange (e.g. Ivanov & Papaloizou 2004;Naoz et al. 2016;Michaely & Perets 2016;Klencki et al. 2017;Wu 2018). In these cases, dissipation in the deep interior would quickly lead to runaway tidal inflation (Bodenheimer et al. 2001;Gu et al. 2003Gu et al. , 2004. By contrast, dissipation exclusively on the surface allows for the possibility of an unperturbed interior of the object, even as large amounts of energy are exchanged (e.g. Wu 2018). Our results are indicate that wave-breaking dissipation occurs exclusively in a surface layer, and further, that this overlying layer damps the oscillations of interior without heat dissipating into the interior.
It is useful to compare these properties to various kinds of linear tidal dissipation. We have so far considered objects with convective envelopes, but when a radiative envelope is present, the principle dissipation mechanism is radiative damping of oscillations. Here radiative losses near the object surface deplete mode energy. These losses are occur because of the low optical depth in these outermost regions (though we note that this dissipation can also occur at an interior radiativeconvective boundary). For low-frequency g-modes, it can be the case that γ rad > ω α , which implies complete damping of the mode every cycle, and a continuous torque as it is re-excited. Zahn (1975) shows that in this limit, a crucial factor is the overlap integral between the mode and the tidal potential, which is generally very small for high-order g-modes. For fundamental modes, while the overlap integral is large, dissipation is inefficient, such that γ rad ω α . One way to understand this result is that the mode eigenfunctions are such that there is significant mode energy in objects' opaque deep interiors, while sufficiently high-order g-modes might concentrate the mode energy closer to the surface (e.g. Fuller & Lai 2012a). Thus, although radiative damping results in surface dissipation of mode energy, it is thought to be an inefficient driver of tidal orbit change in many circumstances.
By contrast, dissipation of oscillatory motions into the random velocity field of convection can be very efficient, especially for modes that couple strongly to the tidal potential, like fundamental modes. In this case, the turbulent viscosity dissipates mode energy at a depth that depends on the mode eigenfunction, the density profile, and the turbulent viscosity as a function of radius (Zahn 1977;Sun et al. 2018;Vick & Lai 2020, equation 19). The turbulent viscosity ν ∼ vH where, v and H describe the characteristic velocity and length scale of eddies. Because the relevant stellar properties of the sound speed and pressure scale height decrease near a stellar surface, convective dissipation is maximized within the object interior, where radiation must diffuse outward over the thermal timescale. The diffusion time of the most dissipative layers and the adjustment of the object remain topics of intense interest (e.g. Podsiadlowski 1996;Barker & Ogilvie 2010;Barker 2011;Ogilvie 2014;Wu 2018;Su et al. 2020). Turbulent dissipation can be very strong, however, with estimated values of γ turb ∼ (M conv /M )(M R 2 /L) −1/3 . Scaled to the dynamical timescale of the object, For example, Vick & Lai (2020) and Vick et al. (2021) report γ turb ∼ 10 −5 t −1 1 for the sun's relatively low-mass mild convective layer (∼ 0.1M ), while γ turb ∼ 0.04t −1 1 for the sun once it ascends to the tip of the giant branch, or for the vigorous convection of near-Eddington massive stars with fully convective envelopes.
We can therefore conclude that for fundamental modes, radiative damping is inefficient relative to nonlinear damping driven by wave breaking, and γ rad γ nl . However, for giants and supergiants with deep convective envelopes, it can be the case that γ turb ∼ γ nl . This comparison carries implications for the occurrence of wave breaking and the evolution of systems under the combination of linear and nonlinear dissipation. We discuss these properties in the following two subsections.

Conditions for Wave Breaking
Our discussion so far has focused exclusively on our model stars, but the lessons derived can be extended to fundamental modes in a variety of stellar and planetary bodies. We find that the criterion |k r ξ r | > 1, evaluated at the object surface at periapse, reliably determines whether wave breaking occurs. By contrast, we note that wave breaking of lower-frequency g-modes can happen in the deep interior of an object (Barker & Ogilvie 2010;Barker 2011), or at the surface (Fuller & Lai 2011, 2012b,c, 2013Vick et al. 2017), depending on an object's structure. In terms of typical stellar parameters, and again adopting k r (R 1 ) ∼ ω/c s (R 1 ), this implies wave breaking for a mode amplitude of where T eff is the surface effective temperature and T vir = (µm p /k b )GM 1 /R 1 is the Virial temperature and, for the sun T eff ≈ 4×10 −4 T vir . In terms of the energy deposited by tides into a fundamental mode, Where, for context, we note that the l = 2 mode masses of our model polytropes are 5 × 10 −4 M 1 and 10 −2 M 1 , respectively for our Γ s = 1.35 and Γ s = 5/3 cases (see Section 3). Thus, depositing a tiny fraction of the object's binding energy into dynamical tidal oscillations can be sufficient to drive wave breaking at the surface of realistic stars. Which systems evolve into wave breaking depends on their orbital and stellar evolution, which are coupled through tidal dissipation (for a related discussion of dissipation and chaos, see Vick & Lai 2018). When γ lin P orb 1, dissipation removes only a small fraction of oscillation energy each orbit and modes can accumulate over many orbital cycles. Here γ lin represents the relevant linear dissipation mechanism. In this way, a star with either a radiative or small convective envelope might experience a mode energy that grows over many orbital cycles until E α ∼ E nl . At that point, wave breaking occurs, and the nonlinear dissipation rate γ nl becomes relevant. Because γ nl P orb 1 for many systems, the mode energy saturates and subsequent orbital evolution is driven by nonlinear wave breaking. By contrast, in systems with vigorous, fully-convective envelopes, γ lin ∼ γ nl . This indicates that mode energies saturate at a level dictated by γ lin P orb and do not necessarily grow over many orbits. This implies that when nonlinear wave breaking occurs in these systems, it is because the single periapse passage energy deposition exceeds the nonlinear wave breaking condition, ∆E α E nl , which is the case in our hydrodynamic simulations. This scenario occurs in systems experiencing rapid stellar radius evolution, such that they remain highly eccentric while evolving toward stronger tidal passages and, eventually, mass exchange (Vick et al. 2021).

Orbital evolution driven by nonlinear wave breaking
As discussed in the previous subsections, wave breaking has the crucial property of dissipating heat at into the atmosphere surrounding an object, rather than its interior. This feature of wave breaking may be crucial in allowing objects to survive tidal circularization with nonlinear dissipation, even when the total energy dissipated is larger than the binding energy of the star or planet. We find that mass loss to the wave breaking layer is an additional important process. As a system evolves through many orbits, we can imagine wave-breaking atmospheres forming when E α E nl . Because the atmosphere extends several times the object's radius, it seems likely that a large fraction will be stripped in the subsequent periapse passage. When nonlinear dissipation is the primary driver of orbital evolution, the rate of change of the mode energy, γ nl , is accompanied by a rate of change of the object mass, through the formation and subsequent removal of wave breaking layers.
The comparison of Figure 9 provides crucial guidance. We observe that for two different polytropic envelope structures, there is an approximately linear relation between E α and m wb , as described by equation (10). In the limit where nonlinear dissipation is the dominant mechanism (γ nl γ lin ) we can then express the relationship between orbital energy change and stellar mass change as which allows one to estimate the fractional mass-loss implied for the orbital energy to change by of order itself, and where C wb ≈ 0.25 in our simulations, equation (10). The relationship between mass loss and orbital change is therefore crucial to the outcome of the nonlinear dissipation phase. For sufficiently small C wb or q, equation (14) shows that the orbital energy can undergo large changes as the object loses small amounts of mass. If C wb or q are larger, we can have |d ln M * /d ln E orb | 1, and the nonlinear phase is dominated by the mass loss from the wave breaking layer. Because many stars and planets become less dense upon mass loss, mass-loss driven evolution can lead to runaway tidal excitation and eccentric Roche lobe overflow rather than circularization. This is precisely the behavior that was observed in the simulations of nonlinear tidal excitation on giant planets by Guillochon et al. (2011). In these models, tides evolved to nonlinear amplitudes, through phases of fundamental mode wave breaking (see their Figure  4, for example), and into a runaway of increasing tidal excitation and mass loss (Guillochon et al. 2011).

Radiative Cooling and Observable Signatures
Wave breaking converts tidal energy to heat through the disordered motions of turbulence. We can estimate the radiative diffusion times through the wavebreaking atmosphere on the basis of their optical depth, τ ≈ ρκL, where L ∼ R 1 is the size scale of the atmosphere, and κ is the opacity. As a simple example, we take M 1 = M , R 1 = R , and κ ≈ 0.34 cm 2 g −1 .
Adopting the approximation, ρ ∼ m atm /R 3 1 . Under these conditions, for a relatively strong encounter with m atm = 10 −4 M 1 , we find τ ∼ 10 7 , and the diffusion timescale, t cool ∼ τ R 1 /c ∼ 1 yr, or approximately 10 4 times the dynamical time of the object. The conclusion that we can draw from this estimate is that for strong encounters, radiative diffusion effectively cools atmosphere layers on timescales that are long compared to the periapse passage and subsequent mode damping, but extremely short compared to the nuclear evolution of a star. However, a weaker encounter, such as one where wave breaking occurs only near the photosphere, where τ ∼ 1, might experience cooling on a much shorter timescale. If cooling causes the atmosphere to reintegrate with the stellar surface layers, it could limit the dissipative efficacy of these weakest wave-breaking cases.
The observable signatures of tidal wave breaking will similarly be varied depending on the particular objects involved along with the mass of the wave breaking layer. Wave breaking shock-heats atmosphere material to temperatures comparable to the object's Virial temperature. Thus post-shock temperatures may be of the order of 10 5 K for a 1M , 100R giant, 10 7 K for a main sequence star, and 10 9 K for a white dwarf. In the presence of an optically thin atmosphere, this hot material would represent itself as an ultraviolet or x-ray excess, with luminosity matching the mode decay rate. Scaling the example of our r p = 1.7R 1 case (Figure 8), the dissipation rate peaks at ∆E/τ diss . For a sun-like star, this implies a dissipative power ofĖ ∼ 10 38 erg s −1 , near the Eddington limit luminosity and enough to drive a surface outflow (Quataert et al. 2016). Scaling to our 100R giant, the dissipative power isĖ ∼ 10 33 erg s −1 . Thus, the relative proportion of the dissipative power can represent very different fractions of the nuclear burning power of these different objects (which might be of the order of 10 33 erg s −1 and 10 36 erg s −1 , respectively). Similarly, while the characteristic wavelength that this radiation emerges at will be dependent upon the optical depth. For example, optically thick dissipation in the presence of a dense atmosphere around a giant star could form dust and emerge as an infrared excess.
More work is needed to clarify these many uncertainties. However, the presence of a disturbed, perhaps shock-heated surface and circumstellar material are likely to be universal -if not unique -traces of surface wave breaking. The case of the massive heartbeat binary system MACHO 80.7443.1718 provides one intriguing example. Here, transient emission lines disappear as the star sweeps through periapse and reappear in the aftermath (Jayasinghe et al. 2021), perhaps indicative of nonlinear tidal wave breaking in action in this system.

SUMMARY AND CONCLUSION
We have performed simulations of close, eccentric tidal passages involving a polytropic model star and a perturber of one tenth its mass. Such a scenario might arise as a consequence of stellar evolution increasing the radius of a star in a binary system, or secular torques in a higher-order multiple system. Indeed, recent work by Vigna-Gómez et al. (2020) and Vick et al. (2021) has emphasized that, especially in systems involving massive stars, the pace of stellar evolution can outstrip that of tidal circularization implying that many such systems approach the onset of mass exchange or a common envelope phase with high eccentricities. Our simulations are relevant to what happens following that phase of gradual evolution into encounters of increasing strength. Some key findings of our study are: 1. In sufficiently close tidal passages, fundamental modes can become large enough in amplitude to lose phase coherence and break, shock heating the outermost layers of the tidally-perturbed star (Figure 3). Wave breaking dissipates the coherent energy of the oscillatory mode into the disordered motion of turbulence (Figure 4), and eventually, heat.
2. The typical nonlinearity condition on the local wave number and amplitude, |k r ξ r | 1, evaluated at periapse, provides a good description of the mass involved in wave breaking ( Figure 5). This mass forms a shock-heated atmosphere layer around the star (Figures 6 and 7).
3. Ongoing oscillations in the stellar interior are damped as they steepen into spiral shock waves in the atmosphere and dissipate ( Figure 6). The nonlinear dissipation rate is set by the rate that these spiral shocks sweep through the atmosphere mass (Figure 8), and thus depends linearly on the atmosphere mass, equation (9).
4. Compared to linear dissipation mechanisms for fundamental modes, nonlinear dissipation by wave breaking is similarly efficient to turbulent dissipation in the fully convective envelopes of giant and supergiant stars, and much more efficient than radiative dissipation. Wave breaking removes mass by lofting it into the atmosphere layer, but has the crucial property of dissipating energy only on the surface, leaving the interior unperturbed.
5. Systems evolving primarily under the influence of nonlinear dissipation lose mass to wave breaking at a rate that is linearly proportional to the total energy change, equation (14). This relationship can predict the outcome of these systems' evolution under the influence of nonlinear dissipation. When |d ln M * /d ln E orb | 1, wave breaking acts to circularize the system orbit during the nonlinear phase (e.g. when q 1). When |d ln M * /d ln E orb | 1, runaway tidal mass loss a likely outcome (e.g. when q 1).
There are many directions for future study. For example, different orbital properties or rotation of the tidal star can change the phase dependence of whether a mode propagates faster or slower than the forcing applied by a perturber. Linear theory predicts different mode spectra depending on system mass ratio, orbital eccentricity, and rotation. Although the general process of oscillation growth up to breaking at a nonlinearity threshold will hold in these systems, the particular modes and their dissipation might differ. Because saturation at the wave-breaking limit seems to be a characteristic feature of the evolution of systems over many orbits, it would be useful to simulate the transition between the breaking and non-breaking tide with more fidelity (e.g. as shown in Figure 9). Radiative cooling of the surface dissipation layers over the course of an orbital period may be important in the dissipative dynamics and fractional mass loss from the wave breaking atmosphere, especially in cases of a thin wave-breaking atmosphere near the nonlinearity limit. Relatedly, there are a broad array of astrophysical systems in which fundamental mode tidal wave breaking may shape the resultant orbital evolution. In particular, systems in which the overall change in orbital energy needed to circularize exceeds the binding energy of the objects themselves -including eccentric formation channels for hot Jupiters and the black hole low mass x-ray binaries -might be particularly sensitive to the details and outcome of nonlinear dissipation.

REPRODUCTION SOFTWARE AND DATA
Software and data to reproduce the results of this study are publicly available in three repositories. The Athena++ simulation setup is available at (MacLeod 2022a), post-processing and analysis software that reproduces the figures is available at (MacLeod 2022b) and selected data needed to reproduce the figures along with simulation runtime parameters are available at (MacLeod 2022c).