\section{Introduction} 
\label{sec:Intro} 

The Advanced Laser 
 Interferometer Gravitational Wave Observatory\footnote{\url{http://www.ligo.caltech.edu/}} (aLIGO) 
 detectors \citep{Aasietal2015} and Advanced Virgo\footnote{\url{http://www.ego-gw.it/}} (AdV) \citep{Acerneseetal2015} 
 have made the first six detections of GWs from approximately circular inspiraling binaries
 \citep{Abbottetal2016c,Abbottetal2016b,Abbottetal2017,LIGOColl2017,Abbottetal2017_2,Abbottetal2017_3}, 
 and opened a new window through which to observe the universe. These advanced gravitational-wave (GW)
 detectors together with upcoming instruments KAGRA\footnote{\url{http://gwcenter.icrr.u-tokyo.ac.jp/en/}} 
 \citep{Somiya2012} and LIGO-India\footnote{\url{http://www.gw-indigo.org/}} \citep{Iyeretal2011,Abbottetal2016a} 
 are expected to continue to detect GW sources in the upcoming years \citep{Abbottetal2016a}. The orbital
 eccentricity was neglected in the analysis of the detected GW sources, but a preliminary upper limit was 
 claimed to be $e \lesssim 0.1$ at 10 Hz \citep{Abbottetal2016e,Abbottetal2016d,Abbottetal2017,Abbottetal2017_4}.
 In this paper we estimate the future potential of the aLIGO-AdV-KAGRA network of advanced GW detectors
 to measure the orbital eccentricity and other physical parameters of initially highly eccentric sources.
 
 Initially highly eccentric black hole (BH) binaries are inspiraling systems, which have orbital eccentricities
 beyond $e_0 \geq 0.9$ when their peak GW frequency \citep{Wen2003} enters the sensitive frequency band of 
 advanced Earth-based GW detectors. The orbital eccentricity decreases in the inspiral phase from this value 
 until the last stable orbit (LSO) \citep{Peters1964}. Such systems can form in multiple ways including 
 single-single encounters due to GW emission \citep{Kocsisetal2006,OLearyetal2009,Gondanetal2017} in 
 dense, high-velocity-dispersion environments; dynamical multibody interactions
 \citep{Gultekinetal2006,OLearyetal2006,Kushniretal2013,AmaroSeoaneChen2016,AntoniniRasio2016,Rodriguezetal2017}; 
 the secular Kozai-Lidov mechanism
 \citep{Wen2003,Thompson2011,Aarseth2012,AntoniniPerets2012,Antogninietal2014,Antoninietal2014,Antoninietal2016,Breiviketal2016,Rodriguezetal2016a,Rodriguezetal2016b,VanLandinghametal2016,Hoangetal2017,PetrovichAntonini2017,SilsbeeTremaine2017,RandallXianyu2018} 
 in hierarchical triples, or the binary-single interaction \citep{Samsingetal2014,SamsingRamirezRuiz2017,Samsing2017,Samsingetal2017}.
 Eccentric BH binaries offer promising new detection candidates.
 
 Previous parameter estimation studies of stellar-mass compact binaries have mostly focused
 on circular binaries 
 (see \citealt{Finn1992,FinnChernoff1993,Markovic1993,CutlerFlanagan1994,JaranowskiKrolak1994,
 Kokkotasetal1994,Krolaketal1995,PoissonWill1995} for the first papers and 
 \citealt{Chatziioannouetal2014,Favata2014,Mandeletal2014,OShaughnessyetal2014,Rodriguezetal2014,Canizaresetal2015,Veitchetal2015,Berryetal2015,Milleretal2015,Farretal2016,Mooreetal2016,Langeetal2017,Vitaleetal2017} for recent developments) 
 due to their predicted high detection rates. The current detections constrain the merger rate 
 density of BH-BH mergers in the Universe to $12-213 \, \mathrm{Gpc}^{-3} \mathrm{yr}^{-1}$ \citep{Abbottetal2017}, 
 which corresponds to a detection rate between $400-7,000 \, \mathrm{yr}^{-1}$ for a typical $2 \, \mathrm{Gpc}$ 
 detection range for aLIGOs design sensitivity. Furthermore, see \citet{Abadieetal2010} for a partial
 list of historical compact binary coalescence rate predictions, and
 \citet{Dominiketal2013,Kinugawaetal2014,Abbottetal2016e,Abbottetal2016h,Abbottetal2016f,Abbottetal2016g,Abbottetal2016j,Belczynskietal2016,Rodriguezetal2016a,Bartosetal2017,McKernanetal2017,Hoangetal2017,Stoneetal2017} and references therein 
 for recent rate estimates.
 
 However, several theoretical studies have shown that the detection rates of highly eccentric BH 
 binaries may be non-negligible. For sources formed by GW-emission in galactic nuclei (GNs), the expected
 aLIGO detection rate at design sensitivity may be higher than $\approx 100 \rm \, yr^{-1}$ if the BH 
 mass function extends to masses above $25 \, \mathrm{M}_{\odot}$ \citep{OLearyetal2009,KocsisLevin2012}. Recently, 
 such heavy BHs have been observed in several LIGO/VIRGO detections 
 \citep{Abbottetal2016h,Abbottetal2017,Abbottetal2017_2}. Additionally, the expected merger rate densities
 in the Kozai-Lidov channel are $1-1.5 \,{\rm Gpc}^{-3} \, {\rm yr}^{-1}$ for BH binaries forming in nuclear
 star clusters without supermassive BHs (SMBHs) through multi-body interactions \citep{AntoniniRasio2016} 
 and $0.14-6.1 \,{\rm Gpc}^{-3} \, {\rm yr}^{-1}$ in isolated triple systems \citep{SilsbeeTremaine2017}. 
 Of order $1-5\,{\rm Gpc}^{-3} \, {\rm yr}^{-1}$ merger rate density is expected for BH 
 binaries forming via the Kozai-Lidov mechanism in globular clusters  \citep{Antoninietal2014,Antoninietal2016,Rodriguezetal2016a}
 and in GNs \citep{AntoniniPerets2012,Hoangetal2017}, and non-spherical nuclear star clusters may 
 produce BH binary merger rates of up to $15\, {\rm Gpc}^{-3} \, {\rm yr}^{-1}$ \citep{PetrovichAntonini2017}.
 Smaller size GNs with intermediate mass BHs may produce higher rates \citep{VanLandinghametal2016}.
 Binary-single gravitational interactions may greatly
 increase the rates \citep{Samsingetal2014,SamsingRamirezRuiz2017,Samsing2017,Samsingetal2017}. 
 In a companion paper \citep{Gondanetal2017}, we have shown that GW capture sources in 
 galactic nuclei, which appear to be circular to within $e<0.2$ near the LSO, may be highly 
 eccentric at the beginning of the detected waveform at $10 \, \mathrm{Hz}$, and that heavier
 BH binaries are expected to be systematically more eccentric in this channel. The ongoing development of detectors 
 towards their design sensitivity at low frequencies may open the possibility of detecting 
 eccentricity in such systems.
 
 In this paper, we determine the expected accuracy with which a network of ground-based interferometric
 GW detectors may determine the physical parameters that describe highly eccentric BH binaries in 
 comparison to circular sources. We investigate how signal-to-noise ratios (SNRs) and parameter 
 measurement errors depend on the initial orbital parameters, particularly the initial pericenter 
 distance and eccentricity.  We examine if it is possible to measure the initial binary parameters 
 (initial eccentricity and pericenter distance) at formation for sources that form in the GW frequency 
 band of the instrument.
 
 Previous GW parameter estimation accuracy studies for eccentric waveforms were carried out for extreme
 mass ratio (EMRI) sources around SMBHs for LISA (\citealt{BarackCutler2004,PorterSesana2010};
 \,\mbox{\citealt{CornishKey2010}};\,\mbox{\citealt{Mikoczietal2012,Nishizawaetal2016}}) and for 
 low-eccentricity stellar-mass compact binaries for Earth-based GW detector network \citep{Sunetal2015}. 
 The premerger localization accuracy of eccentric neutron star binary systems was determined by 
 \citet{KyutokuSeto2014}, and the source localization accuracy was investigated for low-eccentricity 
 binaries by \citet{Maetal2017}.
 
 The parameter space of an eccentric spinning binary waveform is generally very large, 17-dimensional 
 \citep{Vecchio2004,CornishKey2010}. Therefore, state-of-the-art methods such as Monte Carlo Markov 
 Chain calculations (see \citealt{OShaughnessyetal2014} and references therein) are numerically 
 prohibitively expensive to explore the full range of source parameters for a large set of binaries. 
 For Gaussian noise and a large SNR, the posterior distribution function of the measured parameters 
 is generally well approximated by a multidimensional Gaussian, and the parameter measurement errors 
 can be estimated accurately and very efficiently using the Fisher matrix method 
 \citep{FinnChernoff1993,CutlerFlanagan1994,CutlerVallisneri2007}. Using this technique, we determine 
 the physical parameters' measurement accuracy. 

 We restrict this first study to waveforms introduced by \citet{MorenoGarridoetal1994} and 
 \citet{MorenoGarridoetal1995}, which account for part of the leading order post-Newtonian 
 correction, the GR pericenter precession (hereafter simply precession) and neglect other 
 first post-Newtonian and higher order corrections including those due to spins. Future 
 extensions of this work should include higher order post-Newtonian and merger waveforms 
 (see \citealt{Levinetal2011,Csizmadiaetal2012,Eastetal2013} for waveform generators,  \citealt{Damouretal2004,Memmesheimeretal2004,Konigsdorfferetal2005,Konigsdorfferetal2006,Yunesetal2009,TessmerSchafer2010,TessmerSchafer2011,Huertaetal2014,Mikoczietal2015,Mooreetal2016,Tanayetal2016,Boetzeletal2017,CaoandHan2017,HindererBabak2017,Huertaetal2017,Huertaetal2017_2,LoutrelYunes2017} for analytic waveform models, and \citealt{Hinderetal2008,Eastetal2012,Goldetal2012,GoldBrugmann2013,Eastetal2015,Paschalidisetal2015,Eastetal2016,Lewisetal2017} for waveforms of numerical relativity simulations of eccentric compact binary inspirals). In this paper, we focus 
 on BH-BH binaries, but the method is also applicable for neutron star-neutron star and
 neutron star-black hole binaries on highly eccentric orbits as long as tidal interactions
 and matter exchange among the components are negligible (see
 \citealt{Goldetal2012,Eastetal2015,Eastetal2016,Radiceetal2016} and references therein).
 \footnote{Eccentric neutron star (NS) binaries (NS-NS or NS-BH) will also benefit from 
 additional information if an electromagnetic counterpart is identified, which may lead 
 to smaller parameter errors \citep{Radiceetal2016}.}

 Once a large number of GW sources is detected, the correlations between the orbital eccentricity, 
 binary total mass, reduced mass, and spins may be distinctive among different astrophysical mechanisms
 leading to BH mergers \citep{OLearyetal2009,Cholisetal2016,Rodriguezetal2016a,Rodriguezetal2016b,Chatterjeeetal2017,Gondanetal2017,SamsingRamirezRuiz2017,SilsbeeTremaine2017,Kocsisetal2017}. Therefore, detections of eccentric BH binaries have a potential in constraining
 GW source populations. 
 
 However, detecting eccentric sources and recovering their physical parameters is very challenging.
 So far three search methods were developed to find the signals of stellar-mass eccentric BH 
 binaries in data streams of GW detectors \citep{Taietal2014,Coughlinetal2015,Tiwarietal2016}. All
 three methods achieve substantially better sensitivity for eccentric BH binary signals than existing 
 localized burst searches or chirp-like template based search methods. Once a source is detected, 
 different algorithms are used to recover its physical parameters. For compact binary coalescences,
 \textsc{BAYESTAR} \citep{SingerPrice2016} is an online fast sky localization algorithm that produces 
 probability sky maps, \textsc{LALInference} \citep{Veitchetal2015} is an offline full parameter 
 estimation algorithm, and \textsc{gstlal} \citep{Cannonetal2012,Priviteraetal2014} is a low-latency 
 binary BH parameter estimation algorithm. All three algorithms use waveform models of compact binaries 
 on circular orbits. In addition, for short-duration GW ''bursts'' with poorly modeled or unknown waveforms, 
 \textsc{Coherent WaveBurst} \citep{Klimenkoetal2016}, \textsc{BayesWave} \citep{CornishLittenberg2015}, 
 and \textsc{LALInferenceBurst} \citep{Veitchetal2015} pipelines produce reconstructed waveforms with 
 minimal assumptions on the waveform morphology. The development of algorithms recovering the parameters
 of compact binaries on eccentric orbits are currently underway. These algorithms will play an important
 role for the astrophysical interpretation of eccentric sources.
 
 The paper is organized as follows. In Section \ref{sec:WaveMod}, we summarize the basic 
 formulae describing the time-domain and frequency-domain eccentric waveform model. In Section
 \ref{sec:GWdetectors}, we outline the properties of advanced detectors we use in the analysis. 
 In Section \ref{sec:FisherForm}, we describe the signal parameter measurement estimation method. 
 In Section \ref{sec:ParamSpace}, we discuss which parameters of an eccentric binary can be 
 measured through the binary's waveform. We present our main results in Section \ref{sec:NumRes},
 and compare or results with previous papers. Finally, we summarize our conclusions in Section 
 \ref{sec:Conc}. Several details about our methodology is included in appendices. In Appendix 
 \ref{sec:circlim}, we consider the values of source parameters in the circular limit. Next, 
 in Appendix \ref{sec:RoGBD}, we introduce the geometric conventions we use to describe 
 how the GWs interact with ground-based detectors. In Appendix \ref{sec:TimeinDetBand}, we 
 discuss the applicability of the assumptions of neglecting the Earth's rotation around its 
 axis and Earth's motion around the Sun. In Appendix \ref{sec:NumEffSNRFisher}, we derive 
 numerically effective formulae to reduce the computational cost of numerical calculations
 of the $\mathrm{SNR}$ and the Fisher matrix. In Appendix \ref{sec:TestCodes}, we present 
 numerical comparisons to validate our codes for both precessing and non-precessing waveforms.
 
 We use $G=1=c$ units when referring to the initial orbital parameters, and when determining the 
 phases of waveforms. We work in the observer frame assuming a binary at cosmological redshift $z$.
 In this frame, all of the formulae have redshifted mass parameters \mbox{$m_z=(1+z)m$}.\footnote{Additional corrections are necessary if the binary has a peculiar velocity \citep{Kocsisetal2006a}.}
 
 
 
 
\section{Eccentric Waveform Model}
\label{sec:WaveMod}
 
 In this section, we summarize the basic formulae describing the time-domain (Section 
 \ref{sec:TimeWaveMod}) and frequency-domain (Section \ref{subsec:FreqDom}) eccentric waveform models
 including precession in the leading quadrupole-order radiation approximation using the Fourier-Bessel
 decomposition. Note that we neglect the radiation of higher multipole orders, which are typically 
 subdominant at least in cases where the initial pericenter distance is not close to a grazing or 
 zoom-whirl configuration and the initial velocity is much less than the speed of light
 \citep{Davisetal1972,Bertietal2010,Healyetal2016}.
 
 
\subsection{The waveform in time domain} 
\label{sec:TimeWaveMod}
 
 We adopt the waveform model of \citet{MorenoGarridoetal1994} and \citet{MorenoGarridoetal1995}, which 
 describes the quadrupole waveform emitted by a spinless binary on a Keplerian orbit undergoing slow 
 precession. For a fixed semi-major axis $a$ and orbital eccentricity $e$, the two polarization states
 of a GW, $h_+$ and $h_{\times}$, with component masses $m_A$ and $m_B$, and at luminosity distance 
 $D_\mathrm{L}$, can be given in the observer's time-domain as \citep{MorenoGarridoetal1995}: 
\begin{align} \label{eq:hplusBessel}
 \nonumber h_+(t) & = - \frac{ h \, \sin^2 \Theta }{ 2 } \sum_{n=1}^{\infty} A_n \cos \Phi_n(t)
  \\ 
 & + \frac{ h (1+\cos^2 \Theta) }{ 2 } \sum_{n=1}^{\infty} \left( B^+_n \cos \Phi^-_n(t)
  - B^-_n \cos \Phi^+_n(t) \right) \, ,\\
 \label{eq:hxBessel}
  h_{\times} (t) &= - h \, \cos \Theta \sum _{n=1}^{\infty}
  \left( B^-_n \sin \Phi^+_n(t) + B^+_n \sin\Phi^-_n(t)   \right) \, ,
\end{align}
where $\Theta$ is the angle between the orbital plane and the line-of-sight to the 
 observer, $\Phi_n^{\pm}$ describe the orbital phase given below for the $n^{\rm th}$ harmonic,
\begin{equation}  \label{eq:h}
  h = \frac{ 4 \, M_{{\rm tot},z} \, \mu_z }{a \, D_\mathrm{L}  }  \, , 
\end{equation} 
 where $M_{\mathrm{tot},z}= (m_{A}+m_B)(1+z)$ is the redshifted total binary mass at cosmological 
 redshift $z$, and \mbox{$\mu_z =(1+z) m_A m_B (m_A + m_B)^{-1}$} is the redshifted reduced mass. 
 We can express the luminosity distance for a flat $\Lambda$CDM cosmology as a function of $z$ as 
\begin{equation} \label{eq:CovDLz} 
  D_\mathrm{L} = \frac{ (1+z)c }{ H_0 } \int_0 ^z \frac{dz'}{\sqrt{ \Omega_\mathrm{M} 
  (1+z')^3 + \Omega_{\Lambda} }} \, , 
\end{equation} 
 where  $H_0 = 68 \mathrm{km} \, \mathrm{s}^{-1} \mathrm{Mpc}^{-1}$ is the Hubble constant, 
 and $\Omega_\mathrm{M} = 0.304$ and $\Omega_{\Lambda} = 0.696$ are the density parameters 
 for matter and dark energy, respectively \citep{PlanckCollaboration2014a,PlanckCollaboration2014b}. 
 
 The  $A_n$ and $B_n^{\pm}$ prefactors in Equations (\ref{eq:hplusBessel}) and (\ref{eq:hxBessel})
 are the linear combinations of Bessel function of the first kind, $J_n(x)$,
 \begin{equation}\label{eq:AnBn}
  A_n = J_n (ne)\,, \quad B^{\pm}_n = \frac{S_n \pm C_n }{2} \, ,
 \end{equation}
  where $e$ is the orbital eccentricity,
\begin{align} 
  S_n & = - \frac{(1-e^2)^{1/2}}{e} \frac{2}{n} J'_n (ne) + \frac{(1-e^2)^{3/2}}{e^2} 
  2 n J_n (ne) \, , \label{eq:Sn} \\ 
  C_n & = - \frac{ 2-e^2 }{ e^2 } J_n (ne) + \frac{ 2 (1-e^2) }{ e } J'_n (ne) \, , \label{eq:Cn} 
\end{align} 
 and $J'_n(ne)$ is the first derivative of $J_n(ne)$ with respect to $e$, which satisfies
\begin{equation} \label{eq:dJnde} 
  J'_n (ne) = \frac{ n }{ 2 } \left[ J_{n-1} (ne) - J_{n+1} (ne)\right] \, . 
\end{equation} 
 We will also need the second derivative of $J_n(ne)$ with respect to $e$ when calculating the 
 Fisher matrix, thus we introduce $J''_n(ne)$ as
\begin{align} \label{eq:J2na} 
 J''_n &= \frac{n^2}{4} \left[ J_{n-2}(ne)-2J_{n}(ne) + J_{n+2}(ne) \right]  \quad{\rm if~}n\geq 2 
 \,, \\
 J''_1 &=  - \frac{J_{1}(e)}{2}-\frac{1}{4} \left[ J_{1}(e)-J_{3}(e) \right] \, .
\end{align} 
 
 The phase functions $\Phi_n(t)$ and $\Phi_n^{\pm}(t)$ in Equations (\ref{eq:hplusBessel}) and
 (\ref{eq:hxBessel}) are
\begin{align} 
   \Phi_n (t) &= \Phi_c - 2 \pi n \int^{t_c}_{t} \nu(t') dt' \, , \label{eq:Phi} \\
  \Phi^{\pm}_n(t) & = \Phi_n(t) \pm 2\gamma(t) \, , \label{eq:Phipm}
\end{align} 
 where the second term in the right hand side is $n$ times the mean anomaly expressed 
 with the time-integral of the redshifted Keplerian mean orbital frequency $\nu$, $\Phi_c$ is
 the phase extrapolated to coalescence time $t=t_c$, and $\gamma$ is the azimuthal angle of the 
 pericenter relative to the $x$ axis of the coordinate system defined by the orbital plane. 
 The redshifted Keplerian mean orbital frequency may be expressed with the dimensionless 
 pericenter distance 
 \begin{equation}
  \rho_\mathrm{p}=\frac{a(1-e)}{M_{\mathrm{tot},z}}
 \end{equation}
(where $a$ is the semi-major axis in the observer frame) as
\begin{equation}  \label{eq:nuKepler} 
 \nu(e,\rho_\mathrm{p}) = \frac{ (1-e)^{3/2} }{2 \pi \rho_\mathrm{p}^{3/2} M_{\mathrm{tot},z}} \,. 
\end{equation} 
 
 For an inspiraling binary, both the eccentricity and the Keplerian orbital frequency evolve in
 time. Assuming quadrupole radiation and adiabatic evolution of orbital parameters, the equations
 of time evolution of $e$ and $\nu$, as seen at some cosmological redshift, can be given to
 leading order as \citep{Peters1964}
\begin{eqnarray} \label{eq:dedt} 
  \dot{e} &=& - \frac{304}{15} \frac{e \mathcal{M}_z^{5/3} (2 \pi \nu)^{8/3}}{(1-e^2)^{5/2}} 
   \left( 1 + \frac{121}{304} e^2 \right) \, , \\
 \label{eq:dnudt} 
 \dot{ \nu } &=& \frac{48}{5 \pi} \frac{ \mathcal{M}_z^{5/3} (2 \pi \nu)^{11/3}}{(1-e^2)^{7/2}} 
  \left(1 + \frac{73}{24} e^2 + \frac{37}{96} e^4 \right) \, , 
\end{eqnarray} 
 where $\mathcal{M}_z = \mu^{3/5}_z M_{\mathrm{tot},z}^{2/5}$ is the redshifted chirp mass, and 
 the overdot denotes a redshifted time-derivative $\dot{x} \equiv dx/dt$. The fraction of the 
 two equations
\begin{equation}  \label{eq:dnude} 
 \frac{d \nu }{de} = - \frac{18}{19} \frac{\nu}{e} \frac{\left(1 + \frac{73}{24} e^2 + 
 \frac{37}{96} e^4 \right)}{(1-e^2) \left(1 + \frac{121}{304} e^2 \right) } \, ,
\end{equation}
 may be integrated as  \citep{Peters1964,Mikoczietal2012}
\begin{equation} \label{eq:nue}
 \nu(e) = c_0/H(e)
\end{equation} 
 where we define
\begin{equation}  \label{eq:He}
  H(e) = e^{18/19}(1-e^2)^{-3/2} \left(1+\frac{121}{304} e^2 \right)^{\frac{1305}{2299}} \, ,
\end{equation}
 and $c_0$ is an integration constant set by the initial condition $\nu(e_0, \rho_{\mathrm{p}0})
 = \nu_0$ or the conditions at the LSO, \mbox{$\nu(e_{\rm LSO}, \rho_\mathrm{pLSO}) = \nu_{{\mathrm{LSO}}}$}
 (see Equations (\ref{eq:defeLSO}) and (\ref{eq:nuLSO}) below). Equation (\ref{eq:nue}) shows that 
 the product $c_0 = \nu(e) H(e)$ is conserved during the evolution. Similarly, it is straightforward 
 to determine the evolution of the dimensionless pericenter distance
\begin{align}  \label{eq:rhop}
 \rho_\mathrm{p} (e) & = \frac{c_1 e^{12/19}}{ M_{{\rm tot},z}^{2/3} (1 + e)}
 \left(1 + \frac{121}{304} e^2 \right)^{\frac{870}{2299}} 
 \nonumber  \\
 & = \rho_\mathrm{p0} \frac{
 e^{12/19}(1 + e)^{-1}\left[1 + (121/304) e^2 \right]^{\frac{870}{2299}}}{
 e_0^{12/19}(1 + e_0)^{-1}\left[1 + (121/304) e_0^2 \right]^{\frac{870}{2299}}}
\nonumber\\&
 = \rho_\mathrm{pLSO} 
 \frac{ e^{12/19}(1 + e)^{-1}\left[1 + (121/304) e^2 \right]^{\frac{870}{2299}}}{ e_{{\mathrm{LSO}}}^{12/19}
 (1 + e_{{\mathrm{LSO}}})^{-1}\left[1 + (121/304) e_{{\mathrm{LSO}}}^2 \right]^{\frac{870}{2299}} } 
\end{align}
 \citep{Peters1964}, where $c_1=(2\pi c_0)^{-2/3}$ and in the second and third lines we expressed
 the evolution with the initial condition $\rho_\mathrm{p0}$ and $e_0$, or the ``final condition'' 
 at the LSO, which satisfies 
\begin{equation} \label{eq:defeLSO}
 \rho_\mathrm{pLSO} = \rho_\mathrm{p} (e_\mathrm{LSO}) = \frac{6 + 2 e_\mathrm{LSO}}{1 +
 e_\mathrm{LSO}} 
\end{equation}
 in the leading order approximation in the test mass geodesic zero spin limit \citep{Cutleretal1994}.
 This shows that the evolution may be parameterized with the single parameter $e_\mathrm{LSO}$, or the
 two parameters $\rho_{p0}$ and $e_0$. Note that for any $e$, the orbital frequency depends only on the
 single parameter $c_0$, which is set uniquely by $e_{\rm LSO}$ and $M_{\mathrm{tot},z}$ as
 \begin{equation}\label{eq:CC}
 c_0 = \frac{1}{2\pi M_{{\rm tot},z} }\frac{(1-e_\mathrm{LSO})^{3/2}H(e_\mathrm{LSO})}{[\rho_
 \mathrm{pLSO}(e_\mathrm{LSO})]^{3/2}}\,.
 \end{equation}
 
 We restrict our interest to the repeated burst \citep{OLearyetal2009,KocsisLevin2012} and eccentric
 inspiral phases of the waveform model between $0 < e_\mathrm{LSO}\leq e\leq e_0\leq 1$ and neglect 
 the merger and ringdown phases in this analysis. The repeated burst phase starts when the binary is
 formed with initial eccentricity $e_0>0.9$ and initial dimensionless pericenter distance $\rho_{p0}$,
 and the eccentric inspiral phase ends when the binary reaches the LSO with eccentricity $e_\mathrm{LSO}$.
 Note that during the evolution $e$ and $\rho_{p}$ both shrink strictly monotonically in time. 
 
 Let us also note for further use, that the Keplerian redshifted orbital frequency at the end of 
 the assumed eccentric inspiral waveform (i.e. at the LSO) is given by Equations (\ref{eq:nuKepler})
 and (\ref{eq:defeLSO}) as
\begin{equation} \label{eq:nuLSO}
  \nu_\mathrm{LSO} = \nu(e_\mathrm{LSO}) = \frac{ 1 }{2 \pi M_{\mathrm{tot},z} } 
   \left( \frac{ 1 - e_\mathrm{LSO}^2 }{ 6 + 2 e_\mathrm{LSO} }  \right)^{3/2} \,. 
\end{equation}

 Precession leads to a time-dependent $\gamma$ in Equation (\ref{eq:Phipm}). Using the analysis in 
 \citet{Mikoczietal2012}, we adopt pericenter precession from the classical relativistic motion, 
 and assume that the adiabatic evolution of the orbital parameters are governed by Equations 
 (\ref{eq:dedt}) and (\ref{eq:dnudt}). The angle of precession for a single eccentric orbit in
 the test particle geodesic approximation around a Schwarzschild BH is
\begin{equation}
  \Delta \gamma = \frac{ 6 \pi M_{\rm tot} }{ a(1-e^2) } \, .
\end{equation}
 Using an adiabatic approximation, we approximate the redshifted precession rate to be constant 
 during the orbit with
\begin{equation}  \label{eq:dgammadt}
 \dot{ \gamma} \approx \frac{\Delta \gamma}{T} = \frac{3(2 \pi \nu)^{5/3} M_{{\rm tot},z}^{2/3}}
 {1-e^2} \,. 
\end{equation}
 
 The phase functions given by Equations (\ref{eq:Phi}) and (\ref{eq:Phipm}), can be calculated 
 from Equations (\ref{eq:dedt}) and (\ref{eq:nue}) as\footnote{For circular orbits the Fourier 
 phase is conveniently parameterized by $\nu$ \citep{CutlerFlanagan1994}. However, for eccentric
 inspirals, since $\nu(e)$, $\dot{\nu}(e)$, and $\dot{e}(e)$ are given analytically in the PN 
 approximation, the phase is more conveniently parameterized by $e$ 
 \citep{OLearyetal2009,Mikoczietal2012}.} 
\begin{equation}  \label{eq:Phint}
  \Phi_n (t) = \Phi_c + 2 \pi n \int^{e(t)}_0 \frac{\nu(e')}{\dot{e}(\nu(e'),e')} de' 
\end{equation}
 \citep{CutlerFlanagan1994}. The phase functions, which arise due to precession, $\Phi^{+}_n(t)$ 
 and $\Phi^{-}_n(t)$, follow from Equations (\ref{eq:Phipm}), (\ref{eq:nue}), (\ref{eq:dgammadt}) 
 and (\ref{eq:Phint})
\begin{eqnarray} \label{eq:Phinpmt}
  \Phi^{\pm}_n(t) &=& \Phi _n(t) \pm 2 \gamma_c \mp 2 \int^{t_c}_{t} \dot{\gamma} (t') dt' \\ 
  \nonumber
   &=& \Phi _c \pm 2 \gamma_c + \int^{e(t)}_{0} \frac{2\pi n\nu(e')\pm 2\dot{\gamma}(\nu(e'),e')}
   {\dot e(\nu(e'),e')}  de' \, ,
\end{eqnarray}
 where $\gamma_c$ is the angle of periapsis extrapolated to coalescence. Note that $\rho_{p0}$, 
 $t_c$, $\gamma_c$, $\Phi_c$, and $e_{\mathrm{LSO}}$ are free parameters of the waveform. Alternatively, we 
 may use the corresponding initial values $e_0$, $t_0$, $\gamma_0$, $\Phi_0$, and $\rho_{\mathrm{p}0}$.
 
 

\subsection{The waveform in frequency domain }
\label{subsec:FreqDom}
 
 Since the expressions defining the SNR and the Fisher matrix are both given in Fourier space 
 (Section \ref{sec:FisherForm}), we construct the Fourier transforms of the waveform 
 \footnote{ We find that modulations due to Earth's rotation around 
 its axis and Earth's motion around the Sun can be neglected because the signal spends relatively
 short time in the advanced detectors' sensitive frequency band (see Appendix \ref{sec:TimeinDetBand}).} 
\begin{align}
 \label{eq:FThpx} \tilde{h}_{+,\times}(f) & = \int_{-\infty}^{\infty} h_{+,\times}(t) e^{2 \pi i f} dt
                  \, , 
\end{align} 
 where $h_+(t)$ and $h_{\times}(t)$ are given in Equations (\ref{eq:hplusBessel}) and (\ref{eq:hxBessel})
 as an infinite sum over orbital harmonics $n$. In the stationary phase approximation each frequency 
 harmonic splits into a triplet due to precession (see Equation (\ref{eq:dgammadt})) $\mathbf{f} \equiv 
 (f_n, f_n^{\pm})$ \citep{MorenoGarridoetal1995,Mikoczietal2012}, where
\begin{align}
 \label{eq:fn} f_n & = n \nu \, , \\ 
 \label{eq:fnpm} f_n^{\pm} & = n \nu \pm \frac{\dot{\gamma}}{\pi} \, ,
\end{align}  
 and the Fourier transform simplifies to 
\begin{align} 
  \label{eq:SPAp}                     
     \tilde{h}_+ ( \mathbf{f} ) & = - \frac{h_0}{4} \sin^2 \Theta
    \sum _{n=1}^{\infty} A_n \Lambda_n e^{i ( \Psi_n - \pi/4)} 
     \nonumber\\&\quad
   - \frac{h_0}{4} (1 + \cos^2 \Theta)  \sum _{n=1}^{\infty} B^+_n 
            \Lambda_n ^- e^{i ( \Psi_n^- - \pi/4)}    
     \nonumber\\&\quad
   + \frac{h_0}{4} (1 + \cos^2 \Theta) \sum _{n=1}^{\infty} B^-_n \Lambda_n ^+  
   e^{i ( \Psi_n^+ - \pi/4)} \, , \\
 \label{eq:SPAx} \tilde{h}_{\times} ( \mathbf{f}) & = - \frac{h_0}{2} \, \cos\Theta
   \sum _{n=1}^{\infty}  B^-_n \Lambda_n ^+ e^{i ( \Psi_n^+ + \pi/4) } 
    \nonumber\\&\quad    
   - \frac{h_0}{2} \, \cos\Theta \sum _{n=1}^{\infty} B^+_n \Lambda_n ^- 
     e^{i ( \Psi_n^- + \pi/4 )}  \, , 
\end{align}
 where 
\begin{eqnarray} \label{eq:h0}
  h_0 &=& \frac{ 4 \mathcal{M}_z^{5/3} (2 \pi \nu)^{2/3} }{ D_\mathrm{L} } \, ,\\\label{eq:Lambda}
    \Lambda_ n  &=& \frac{1}{\sqrt{\arrowvert n \dot{\nu} \arrowvert }}\, ,\quad 
  \Lambda_ n^{\pm} = \frac{1}{ \sqrt{\arrowvert n \dot{\nu} \pm \ddot{\gamma}/ \pi 
  \arrowvert } } \, ,
\end{eqnarray}
 $\dot{\nu}$ and $\ddot{\gamma}$ are given by Equations (\ref{eq:dnudt}), (\ref{eq:dnude}) and 
 (\ref{eq:dgammadt}), and \mbox{$\mathcal{M}_z = (1+z) \mathcal{M}$}. The $ \Psi_n $ and $ \Psi_n^{\pm}$
 phases and their first ($\dot{\Psi}_n\, ,\dot{\Psi}^{\pm}_n$) and second ($\ddot{\Psi}_n\, , 
 \ddot{\Psi }^{\pm}_n$) derivatives with respect to redshifted time $t$ are \citep{Mikoczietal2012}
\begin{align}
  \Psi_n (e, f_n) & = 2 \pi f_n t_n - \Phi_n \, , \label{eq:PsiN}  \\
  \Psi^{\pm}_n (e, f^{\pm}_n) & = 2 \pi f_n^{\pm} t_n^{\pm} - \Phi_n^{\pm} \, .   \label{eq:Psipm} 
\end{align}  
 Here the $(t_n,t_n^{\pm})$ parameters of the stationary phase approximation specify the 
 times at which the orbital frequency satisfies Equations (\ref{eq:fn}) and (\ref{eq:fnpm}) 
 for given $(f_n,f_n^{\pm})$, see Appendices A and B in \citet{Mikoczietal2012} for details.
 In Equations (\ref{eq:PsiN}) and (\ref{eq:Psipm}), $(\Phi_n, \Phi_n^{\pm})$ are to be 
 substituted from Equations (\ref{eq:Phint}) and (\ref{eq:Phinpmt}). We eliminate $\nu(t)$ 
 and $e(t)$ for $(f_n,f_n^{\pm})$ using Equations (\ref{eq:nue}), (\ref{eq:CC}), and 
 (\ref{eq:dgammadt}) together with Equations (\ref{eq:fn}) and (\ref{eq:fnpm}) to obtain 
 the frequency-domain waveform.\footnote{In practice, there 
 are closed analytic expressions for the $e$-dependence of $\nu$, $\dot{\nu}$, $\dot\gamma$, 
 $\ddot\gamma$, $\Phi_{n}$, $\Phi_n^{\pm}$, and hence also for $f_n$ and $f_n^{\pm}$. We 
 must invert these relations $f_n(e)$ and $f_n^{\pm}(e)$ to obtain the waveform in frequency
 domain.} The result depends on constant parameters $t_c$, $\Phi_c$, $\gamma_c$, $e_\mathrm{LSO}$,
 $\mathcal{M}_z$, and $\mathrm{M}_{\mathrm{tot},z}$. Further, we note that if the precessing 
 eccentric BH binary forms with $\rho_{\mathrm{p}0}$ and $e_0$, then the frequency-domain 
 waveform is truncated at the corresponding minimum frequency \mbox{$(f_{n,\min}$, $f_{n,\min}
 ^{\pm})=(n\nu_0,n\nu_0\pm \dot\gamma_0/\pi)$}. Furthermore, the waveform model becomes invalid after reaching the LSO (with $\rho
 _{\mathrm{p LSO}}$ and $e_\mathrm{LSO}$), which corresponds to a maximum frequency for each 
 harmonic \mbox{$(f_{n,\max}$, $f_{n,\max}^{\pm})=(n\nu_{\mathrm{LSO}},n\nu_{\mathrm{LSO}}\pm \dot\gamma_{\mathrm{LSO}}/\pi)$} where this model is applicable. If we 
 truncate the waveform at these maximum frequencies, this respectively introduces an explicit
 $(\rho_{\mathrm{p}0}, e_0)$ and $(\rho_{\mathrm{p LSO}}, e_\mathrm{LSO})$ parameter dependence
 in the waveform model. This is shown in Equation~(\ref{eq:hPPC}) in Appendix~\ref{sec:NumEffSNRFisher}.
 Examples of the frequency-domain waveforms are shown in \citet{{KocsisLevin2012}}.
 
 In principle, the number of spectral harmonics of an eccentric binary system is infinite. Note
 however that a large fraction of the signal power is accumulated in a finite number of harmonics.
 Therefore, in order to reduce the necessary computation time, we truncate $n$ at $n_\mathrm{max}
 (e_0)$ \citep{OLearyetal2009,Mikoczietal2012}
\begin{equation}  \label{eq:nmax}
  n_\mathrm{max}(e_0) = \left\{ 5 \frac{ (1+e_0)^{1/2} }{ (1-e_0)^{3/2} } 
  \right\} \, ,
\end{equation}
 which accounts for $99\%$ of the signal power \citep{Turner1977}. Here the bracket $\{ \}$ 
 denotes the floor function. In Appendix \ref{sec:NumEffSNRFisher}, we discuss other technical
 details to optimize the calculation of the SNR and the Fisher matrix.
 
 To test our calculations, we examine the limiting cases of no precession  ($\dot{\gamma}\rightarrow
 0$) and circular orbits ($e\rightarrow 0$), respectively. In Appendix \ref{sec:NumEffSNRFisher},
 we discuss numerical and analytic tricks to optimize the calculation and discuss results for the
 precessing (\emph{Prec}) and precession-free (\emph{NoPrec}) waveform model (i.e. $\dot{\gamma} 
 \equiv 0$).
  
 
 

\section{GW Detectors Used in the Analysis}
\label{sec:GWdetectors}
 
 Here we summarize the GW detectors and the assumed properties of the detector noise in our analysis.
 
 The aLIGO and AdV detectors completed their first two observing run, and made the first six detections of
 GWs \citep{Abbottetal2016c,Abbottetal2016b,Abbottetal2017,Abbottetal2017_2,Abbottetal2017_3,LIGOColl2017}.
 Two additional GW detectors are planned to join the network of aLIGO and AdV; (i) the Japanese KAGRA is 
 under construction with baseline operations beginning in $2018$ \citep{Somiya2012}; while (ii) the proposed
 LIGO-India is expected to become operational in $2022$ \citep{Iyeretal2011,Abbottetal2016a}. LIGO-India was
 approved by the government of India and a study has already suggested site location and orientations of arms 
 for the detector based on scientific figures of merit \citep{Raffaietal2013}. These parameters, however, have
 not been finalized yet, and because of this we omit LIGO-India from the analysis. 

\begin{table}
\centering  
   \begin{tabular}{@{}lrrr}
     \hline
       Detector  &  East Long.  &  North Lat. &  Orientation $\psi$  \\
     \hline\hline
       LIGO H  & $ -119.4^{\circ}$ & $46.5^{\circ}$ & $-36^{\circ}$ \\
       LIGO L  &   $ -90.8^{\circ}$   &   $ 30.6^{\circ} $ & $-108^{\circ}$  \\
       VIRGO &    $ 10.5^{\circ}$  &  $ 43.6^{\circ}$ & $ 20^{\circ}$  \\
       KAGRA &    $ 137.3^{\circ}$ &   $ 36.4^{\circ}$ & $ 65^{\circ}$  \\
     \hline
   \end{tabular} 
   \caption{ \label{tab:DetCord} Locations and orientations of considered GW detectors in 
   the coordinate system defined in Appendix \ref{sec:RoGBD}. LIGO H marks the Advanced LIGO
   detector in Hanford, and LIGO L marks the Advanced LIGO detector in Livingston. }
\end{table}

 Due to the expected similarities of design sensitivities of the aLIGO, AdV, and KAGRA 
 detectors  within the frequency range of BH inspiral waveforms, for simplicity we adopt
 the design sensitivities of the two aLIGO \citep{Abbottetal2016a} for AdV \citep{Abbottetal2016a}
 and for KAGRA \citep{Somiya2012} detectors. Table \ref{tab:DetCord} gives the locations
 and orientations of these detectors, which we used to calculate the response functions.
 For each detector, we define the detector's orientation angle, $\psi$, as the angle 
 measured clockwise from North between the $x$-arm of the detector (see Appendix 
 \ref{sec:RoGBD} for the geometric conventions of detectors) and the meridian that 
 passes through the position of the detector.
 
 We assume that the noise in each detector is stationary colored Gaussian with zero mean, and 
 that it is uncorrelated between different detectors. In reality, detector noise arises 
 from a combination of instrumental, environmental, and anthropomorphic sources that are difficult 
 to characterize precisely \citep{Aasietal2012,Asoetal2013,Aasietal2015}, and non-Gaussian 
 noise transients (glitches) may arise as well \citep{Blackburnetal2008}. However, there 
 are existing techniques to identify and remove glitches from GW strain channels 
 and to reduce the level of these artifacts 
 \citep{LittenbergCornish2010,Prestegardetal2012,Biswasetal2013,Powelletal2015,Boseetal2016,Torresetal2016,Georgeetal2017,Mukundetal2017,Powelletal2017,Shenetal2017}. 
 Furthermore, correlated noise between widely separated detectors can arise from so-called 
 Schumann resonances (predicted in \citet{Schumann1952a,Schumann1952b}  and observed soon 
 thereafter \citep{SchumannKonig1954,BalserWagner1960}), as well as from other EM phenomena 
 such as solar storms, currents in the van Allen belt \citep{Rycroft2006}, and anthropogenic 
 emission \citep[see][and references therein]{Shvetsetal2010,Thraneetal2013}. Note however, that 
 Schumann resonances mostly affect the stochastic GW background searches \citep{Thraneetal2013},
 and a strategy against such noise artifact already exists \citep{Thraneetal2014}. Our 
 simplifying assumptions on uncorrelated Gaussian noise are therefore partly justified.
 

 
\section{Overview of the Fisher Matrix Formalism} 
\label{sec:FisherForm} 
 
 In this section, we provide a brief overview of the Fisher matrix method to estimate the measurement 
 errors of physical parameters characterizing a precessing eccentric BH binary source, and refer 
 the reader to \citet{Finn1992} and \citet{CutlerFlanagan1994} for further details. 
 
 The output of a GW detector, $s(t)$, is a combination of a signal, $h(t)$, and a noise term, $n(t)$; i.e.
\begin{equation}  \label{eq:st} 
  s(t) = h(t) + n(t) \, .
\end{equation} 
 We assume the noise of a detector to be stationary, zero mean, and Gaussian, where the different
 Fourier components of the noise are uncorrelated, i.e. 
 \begin{equation}  \label{eq:defSa} 
  \left\langle \tilde{n}(f) \tilde{n}^*(f') \right\rangle = \frac{1}{2} \delta(f-f') S_n(f) \,  
\end{equation}
 \citep{Nissankeetal2010}, where $\langle\cdot\rangle$ denotes the average, $S_n(f)$ is the one sided noise 
 power spectral density of the detector, the ${}^*$ superscript denotes complex conjugate. With these assumptions, 
 the probability for the noise to have some realization $n_0(t)$ is given as
\begin{equation}  \label{eq:noiseprob} 
   p\left(n \equiv n_0 \right) \propto e^{- \left( n_0 \mid n_0 \right)/2 } \, 
\end{equation} 
 \citep{Finn1992}, where $p(n)$ is the probability distribution function of the noise to assume
 a value $n$, and $(\ldots | \ldots)$ denotes the following inner product between any two functions
 of frequency, e.g. x(f) and y(f):
\begin{equation}  \label{eq:InnProd} 
  \left( x \mid y \right) \equiv 4 \int_{0}^{\infty} 
   \frac{ \tilde{x}(f) \, \tilde{y}^*(f) }{S_n(f)} \, df \, .
\end{equation} 
 The optimal SNR is given by the standard expression 
\begin{equation}  \label{eq:SNR2_1det} 
  \frac{S}{N} = \sqrt{ \left(h \mid h \right) } = \sqrt{ 4 \int_{0}^{\infty} \frac{|\tilde{h}(f)|^2 
  }{S_n(f)} df } \, . 
\end{equation} 

 Here the signal waveform, $h(f)$, depends on the parameter set $\left\{ \lambda_p \,|\, p \in\{ 1 
 \ldots P \right\} \}$, which characterizes the source. For a large SNR, the parameter estimation
 errors \mbox{$\bm{\Delta \lambda}=\left\{ \Delta \lambda_p \,|\, p \in \{1, \ldots, P \right\} \}$} 
 defined as the measured value minus the true value have the Gaussian probability distribution 
 for a given signal
\begin{equation}  \label{eq:probdistpar}
  p ( \bm{\Delta \lambda} ) = \mathcal{N} \, \exp \left( - \frac{1}{2} 
  \Gamma_{ij} \Delta \lambda_i \Delta \lambda_j \right)
\end{equation}
 \citep{Finn1992}, where $\mathcal{N}$ is a normalization constant. In Equation (\ref{eq:probdistpar}), 
 we assume summation over repeated indices, and $\Gamma_{ij}$ is the Fisher information matrix defined as
\begin{equation}  \label{eq:Fisher}
  \Gamma_{ij} \equiv \left( \partial_i h \mid \partial_j h \right) = 
  4 \int_0^{\infty} \frac{ \Re \left( \partial_i \tilde{h}^*(f) \partial_j 
  \tilde{h}(f) \right) }{ S_n(f) } \, df \, ,
\end{equation}
 where  $\partial_i h = \partial h / \partial \lambda_i$ and $\Re$ labels the real part.
 
 Following \citet{CutlerFlanagan1994}, we define the combined signal-to-noise ratio of 
 a network of detectors ($\mathrm{SNR}_\mathrm{tot}$) as an uncorrelated superposition 
 of individual SNRs 
\begin{equation}  \label{eq:SNRtot}
   \left(\frac{S}{N}\right)^2_{\mathrm{tot}} = \sum_{k=1}^{N_\mathrm{det}} \left(\frac{S}{N}\right)^2_k \, ,
\end{equation} 
 where the number of detectors in the network is denoted by $N_\mathrm{det}$, and $(S/N)_k$
 denotes the SNR in the $k^{\mathrm{th}}$ detector. 

 Similarly, for uncorrelated Gaussian noise, the Fisher matrix of a network of detectors is the
 sum of the Fisher matrices of individual detectors,
\begin{equation} \label{eq:Fishertot}
  \Gamma_{ij,\mathrm{tot}} = \sum _{k=1}^{N_\mathrm{det}} \Gamma_{ij,k} \, .
\end{equation}
 The covariance matrix is defined with the inverse of the Fisher matrix:
\begin{equation} \label{eq:CovMatrix}
  \Sigma_{ij} = \left( \Gamma_{ij,\mathrm{tot}} \right)^{-1} = 
   \langle \Delta \lambda_i \Delta \lambda_j \rangle \, ,
\end{equation}
 where the angle brackets denote an average over the probability distribution function in 
 Equation (\ref{eq:probdistpar}). The root-mean-square parameter measurement error $\sigma_i$
 in the parameters $\lambda_i$ marginalized over all other parameters is
\begin{equation}  \label{eq:Sigma}
  \sigma_i = \langle \left(\Delta\lambda_i \right)^2 \rangle^{1/2} = \sqrt{\Sigma_{ii}} \, .
\end{equation}
 The off-diagonal elements of $\Sigma_{ij}$ give the cross-correlation coefficients 
 between parameters $\lambda_i$ and $\lambda_j$.
 
 Parameters can be measured independently if the corresponding Fisher matrix $\Gamma_{ij,\rm
 tot}$ is nonsingular. Otherwise, if the Fisher matrix is singular, then the eigenvector(s) 
 corresponding to the zero-eigenvalue(s) of the Fisher matrix represent the linear combination(s)
 of the parameters, which cannot be measured by the network.
 
 We derive efficient formulae to compute the SNR and the Fisher matrix in Appendix
 \ref{sec:NumEffSNRFisher}. 
 
 
 
 
\section{Measuring the Parameters of Precessing Eccentric Black Hole Binaries}
\label{sec:ParamSpace} 
 
 In this section, we identify the parameters of a precessing eccentric binary that can be extracted
 from the detected waveform for the signal model introduced in Section \ref{sec:WaveMod}. We set the
 parameters in our calculations and measure their errors as follows.
 \begin{itemize}
  \item $D_\mathrm{L}:$ 
  We set $D_\mathrm{L}=100 \,\mathrm{Mpc}$, and measure its relative error $\langle \Delta D_{\rm L}^2 
  \rangle^{1/2} / D_{\rm L} = \langle (\Delta \ln D_{\rm L})^2 \rangle^{1/2}$. This choice is arbitrary, 
  smaller than the nearest circular BH-BH merger detection to date, $340 \pm 140 \, \mathrm{Mpc}$ 
  \citep{LIGOColl2017}. The Fisher matrix method gives accurate results for the parameter measurement 
  errors for high SNR. For moderately larger distances, the errors scale as $\propto D_{\rm L}$. 
  
  \item $\theta_N$ and $\phi_N$: 
  We generate an isotropic random sample of the sky position angles $\theta_N$ and $\phi_N$ by drawing
  $\cos \theta_N$ and $\phi_N$ from a uniform distribution between $[-1,1]$ and $[0,2 \pi]$, and calculate
  the parameter estimation covariance for each sample. The errors of the sky position is described by 
  a localization ellipse. We characterize the sky localization accuracy either by the corresponding 
  proper angular length of the semi-major and semi-minor axes of the sky-localization error ellipsoid
  given by \citet{LangHughes2006}, $(a_N,b_N)$, or its proper solid angle $\Delta\Omega_N = \sqrt{\pi 
  a_N b_N}$. The calculated results are valid if $a_N\ll 1$ radian and $b_N\ll 1$ radian. 

 \item $\theta_L$ and $\phi_L$: We draw the angular momentum vector direction angles from an isotropic
 distribution and construct their error ellipsoids or solid angles similar to that given for $\theta_N$ 
 and $\phi_N$.

 \item $m_A$ and $m_B$: We fix the fiducial component masses to $m_A=m_B=30 \, \mathrm{M}_{\odot}$, consistent
  with the first discovered source GW150914 \citep{Abbottetal2016b}. Such high mass sources are
  expected in galactic nuclei since mass segregation helps to increase their numbers relative to 
  the lower mass binaries, and the SNR is also higher for these binaries \citep{OLearyetal2009}. 
  Since we neglect additional post-Newtonian
  corrections of the GW phase, we restrict the measurement error estimation to $\mathcal{M}_z$ 
  for calculations evaluated for comparison in which we neglect precession.
  However, generally the assumed precessing eccentric waveform model depends on two independent combinations of component 
  masses: $\mathcal{M}_z$ sets the inspiral rate, and $M_{\mathrm{tot},z}$ sets both the apsidal
  precession rate and the final frequency at the LSO. We calculate the relative errors for both
  of these mass parameters and for the precessing eccentric waveform model \mbox{$\langle \Delta 
  \mathcal{M}_z^2 \rangle^{1/2}/ \mathcal{M}_z = \langle( \Delta \ln \mathcal{M}_z)^2 \rangle^
  {1/2}$}, and similarly for $M_{\mathrm{tot},z}$. 

 \item $t_c$, $\Phi_c$, and $\gamma_0$: These parameters only enter in the complex phase of the
  waveform through $\Psi_n$ and $\Psi_n^{\pm}$ (see Equations (\ref{eq:PsiN}) and (\ref{eq:Psipm})), 
  but do not affect the SNR. Since these parameters are responsible for an overall phase shift of
  the waveform, we do not randomize their values but assume the fiducial value $t_c = \Phi_c = 
  \gamma_c = 0$ for each binary in the Monte Carlo sample.

 \item $e_\mathrm{LSO}$: The adopted eccentric inspiral waveform model depends explicitly on 
  the final eccentricity at the LSO, see Equation (\ref{eq:defeLSO}). This quantity parameterizes 
  the evolutionary path of the binary during its eccentric inspiral in the $(\rho_\mathrm{p},e)$
  plane as shown in Equations (\ref{eq:rhop}) and (\ref{eq:defeLSO}); see also Figure 3 in 
  \citet{KocsisLevin2012} for illustration. In fact, any segment of the evolutionary path 
  $\rho_\mathrm{p}(e)$ specifies the value of $e_\mathrm{LSO}$ uniquely. Conversely, $e_
  \mathrm{LSO}$ specifies $\rho_\mathrm{p}(e)$, which sets a constraint on the possible values
  of $(\rho_{\mathrm{p}0},e_0)$, if the post-Newtonian binary inspiral model is extrapolated 
  backwards in time. Indeed, in some cases this is the only indirect information we may have 
  on the formation parameters $(\rho_{\mathrm{p}0},e_0)$. In particular, $e_0\leq 1$ puts an 
  upper bound on $\rho_{\mathrm{p}0}$ for a given $e_\mathrm{LSO}$.\footnote{When studying 
  the measurement errors for non-precessing eccentric binaries, the waveform depends 
  explicitly on a single combination of $e_\mathrm{LSO}$ and $M_{\mathrm{tot},z}$ parameters
  $c_0$ (Section \ref{sec:TimeWaveMod}). Therefore, we use $c_0$ for the \emph{NoPrec} model 
  to avoid a singularity of the Fisher matrix.}
 
 \item $e_0$: We choose several $e_0$ values from the highly eccentric ($e_0 \geq 0.9$) limit when
  discussing the $e_0$ dependence of measurement errors (see Section \ref{subsec:ResParamDist}). 
  However, we restrict to $e_0 = 0.9$ for calculations of a large survey of binaries.\footnote{We 
  note that the $e_0$ dependence of the waveform is due to the truncation of the time-domain 
  waveform for times when $e<e_0$.}
 
 \item $\rho_{\mathrm{p}0}$:  We examine two values for the dimensionless initial pericenter 
  distance $\rho_{\mathrm{p}0} = \lbrace 10,20 \rbrace$, and the circular limit corresponds to 
  $\rho_{\mathrm{p}0} \rightarrow \infty$ \citep{OLearyetal2009}. These values are likely for 
  sources that form through the GW capture mechanism in high velocity dispersion environments 
  such as GN as shown in \citet{OLearyetal2009,Gondanetal2017} or the core collapsed regions of 
  star clusters without a central massive black hole \citep{Kocsisetal2006,AntoniniRasio2016}.
 \end{itemize} 

 If the peak GW frequency of the initial orbit is large enough to be in the detectors' sensitive 
 frequency band then $\rho_{\mathrm{p}0}$ and $e_0$ are directly measurable due to the truncation
 of the time-domain waveform for times when $e<e_0$ and $\rho_{\mathrm{p}}>\rho_{\mathrm{p}0}$.
 In the opposite case only a lower limit may be given for $\rho_{\mathrm{p}0}$, which corresponds
 to $e_0 \rightarrow 1$ \citep{KocsisLevin2012}.

 In summary, we use the following free parameters in the Fisher matrix analysis:
\begin{align}
   \bm{\lambda}_\mathrm{Prec}=& 
   \lbrace \mathrm{ln} (D_\mathrm{L}), \mathrm{ln} 
  (\mathcal{M}_z), \mathrm{ln} (M_{{\rm tot},z}), \theta_N,\phi_N,\theta_L,\phi_L,e_0, 
  \nonumber\\ & \quad 
  e_\mathrm{LSO}, t_c,\Phi_c,\gamma_c \rbrace \, .
\end{align}
 Given these parameters, other parameters' marginalized measurement errors may be determined
 by linear combinations of the covariance matrix based on Equation (\ref{eq:Sigma}). For example, 
 $\rho_{\mathrm{p}0}$ is given by $e_0$ and $e_\mathrm{LSO}$ using Equations (\ref{eq:rhop}) and
 (\ref{eq:defeLSO}). Its measurement error is 
\begin{align}  \label{eq:deltarhop0}
 \left\langle(\Delta \rho_{\mathrm{p}0})^2 \right\rangle
 = & \left( \frac{\partial \rho_\mathrm{p}(e_0,e_\mathrm{LSO})}{\partial e_0} \right)^2 \left\langle (\Delta e_0)^2
 \right\rangle
   \nonumber\\ 
   & + \left( \frac{\partial \rho_\mathrm{p}(e_0,e_\mathrm{LSO})}{\partial e_\mathrm{LSO}} \right)^2 
  \left\langle (\Delta e_\mathrm{LSO})^2 \right\rangle
  \nonumber\\
  & + 2 \frac{\partial \rho_\mathrm{p} (e_0,e_\mathrm{LSO})}{\partial e_0} \frac{\partial \rho_\mathrm{p}(e_0,
  e_\mathrm{LSO})}{\partial e_\mathrm{LSO}} \left\langle \Delta e_0 \Delta e_\mathrm{LSO} \right\rangle\,.
\end{align}
 The parameter estimation errors of individual component masses or the mass ratio can be estimated 
 similarly using $\Delta M_{\rm tot}$ and $\Delta \mathcal{M}$ after inverting $M_{\rm tot}
 (m_a, m_b)$ and $\mathcal{M}(m_a, m_b)$.



\section{Results}
\label{sec:NumRes}
 
 The measurement errors depend on the sky position of the source with respect to the detectors 
 and on the relative orientation of the binary. We generate random Monte Carlo samples of $\sim 
 4500$ binaries by drawing from isotropic distributions of the sky position and of the binary
 orientation normal vector. We present the results for the $\mathrm{SNR}_\mathrm{tot}$ for 
 detecting precessing highly eccentric BH binaries with the GW detector network described 
 in Table \ref{tab:DetCord} and for the expected parameter measurement errors.
 
 \begin{figure} 
    \centering 
    \includegraphics[width=79mm]{SNR_Prec} 
\caption{  Smoothed probability density function of the total network signal-to-noise ratio 
         ($\mathrm{SNR}_\mathrm{tot}$) of gravitational-wave detection from  $30 \, \mathrm{M}_{\odot}-30 \, 
         \mathrm{M}_{\odot}$ precessing eccentric BH binaries with initial eccentricity $e_0=0.9$ and dimensionless 
         pericenter distance $\rho_{\mathrm{p}0} = 10$ and $20$ (green dashed and red dash-dot) and 
         similar binaries in the circular limit (blue solid) at luminosity distance $D_\mathrm{L} = 
         100 \, \mathrm{Mpc}$ with a random source direction and orientation. Distributions correspond 
         to a Monte Carlo sample of $4500$ binaries. Parameters of the assumed detector network are 
         given in Table \ref{tab:DetCord}. The medians of $\mathrm{SNR}_\mathrm{tot}$ distributions 
         are $108.7$, $202.7$, and $137.3$ in the circular limit and for binaries with $\rho_{\mathrm{p}0} 
         = 10$ and $20$, respectively. Systematically higher $\mathrm{SNR}_\mathrm{tot}$ values for 
         precessing highly eccentric BH binaries implies that they are detectable to a larger distance 
         compared to precessing eccentric BH binaries in the circular limit. } \label{fig:SNR_Prec} 
\end{figure} 
 
 
 
\subsection{Signal-to-noise ratio distributions} 

 Figure \ref{fig:SNR_Prec} displays the distribution of the $\mathrm{SNR}_\mathrm{tot}$ for precessing
 highly eccentric BH binaries detected with the detector network described in Table \ref{tab:DetCord}, 
 assuming binary parameters $m_A=m_B=30 \, \mathrm{M}_{\odot}$, $e_0=0.9$, $\rho_{\mathrm{p}0} = \lbrace 10, 20 
 \rbrace$, and similar binaries in the circular limit (see Appendix \ref{sec:circlim} for details).
 Generally, similar to the results of \citet{OLearyetal2009} (see Figure 11 therein, which corresponds
 to a single aLIGO detector), the $\mathrm{SNR}_\mathrm{tot}$ is systematically higher for binaries with
 $\rho_{\mathrm{p}0} = 10$ than for binaries with $\rho_{\mathrm{p}0} = 20$. We find that increasing 
 the initial eccentricity from $e_0=0.9$ to $0.97$ for fixed $\rho_{\mathrm{p}0}$ does not change the 
 $\mathrm{SNR}_\mathrm{tot}$ significantly (see Table \ref{tab:PrecBine0rho0Dep} below), hence we expect
 that the distribution of $\mathrm{SNR}_\mathrm{tot}$ for fixed $\rho_{\mathrm{p}0}$ converges in the $e_0 
 \rightarrow 1$ limit. This is expected since Figure 10 in \citet{OLearyetal2009} shows that a low amount
 of the $\mathrm{SNR}_\mathrm{tot}$ accumulates near $e_0 \approx 1$ for low to moderately high BH masses.  

\begin{figure}
    \center
    \includegraphics[width=79mm]{SNRtot_rho}
\caption{ The $\mathrm{SNR}_\mathrm{tot}$ of precessing eccentric BH binaries as a function of 
          their initial dimensionless pericenter distance $\rho_{\mathrm{p}0}$. Parameters of
          the assumed detector network are given in Table \ref{tab:DetCord}. Here the luminosity 
          distance is $D_\mathrm{L}=100 \, \mathrm{Mpc}$, the initial eccentricity is $e_0 = 0.9$, 
          the source direction and orientation angular parameters are fixed at $\theta_N=\pi/2$, 
          $\phi_N=\pi/3$, $\theta_L=\pi/4$, and $\phi_L =\pi/5$. Depending on the binary mass, 
          the $\mathrm{SNR}_ \mathrm{tot}$ has a maximum at $\rho_{\mathrm{p}0}$ between $9$ and 
          $20$, and converges asymptotically to the value of precessing eccentric binaries in 
          the circular limit for high $\rho_{\mathrm{p}0}$. Note that we find similar trends 
          with $\rho_{\mathrm{p}0}$ for other random choices of binary direction and orientation 
          (not shown). }  \label{fig:SNRtot_rho}
\end{figure} 

 Figure \ref{fig:SNRtot_rho} shows that the $\mathrm{SNR}_\mathrm{tot}$ increases rapidly with $\rho_{
 \mathrm{p}0}$ for low $\rho_{\mathrm{p}0}$, has a maximum between $\rho_{\mathrm{p}0}\sim 9$ and 
 $20$, and decreases for higher $\rho_{\mathrm{p}0}$ approaching the circular binary limit for $\rho
 _{\mathrm{p}0} \rightarrow \infty$. These findings may be understood qualitatively as follows. 
 Within $\rho_{\mathrm{p}0}\leq 40.9 (M_{\rm tot,z}/20 \, \mathrm{M}_{\odot})^{-2/3}$ the binary forms with a 
 characteristic frequency above $10 \, \mathrm{Hz}$ in the detector band \citep[see Equation 
 (59) in][]{Gondanetal2017}. The rapid decrease of the $\mathrm{SNR}_\mathrm{tot}$ for decreasing 
 $\rho_{\mathrm{p}0}<9$ is due to the fact that we neglect the GWs of the first hyperbolic encounter 
 \citep{Kocsisetal2006}. The decrease of the $\mathrm{SNR}_\mathrm{tot}$ at high $\rho_{\rm p0}$ is 
 due to the fact that part of the GW spectrum falls outside of the detectors' sensitive frequency 
 band. For very large $\rho_{\mathrm{p},0}$, the binary becomes circular by the time it enters the 
 detectors' sensitive frequency band, and a significant fraction of the $\mathrm{SNR}_\mathrm{tot}$ 
 accumulates only in the $n=2$ harmonic (Figure \ref{fig:time}), which explains the flat asymptotics 
 for high $\rho_{\mathrm{p}0}$. Decreasing $\rho_{\mathrm{p}0}$ from high to moderate values higher 
 harmonics start to contribute to the $\mathrm{SNR}_\mathrm{tot}$ (Figure \ref{fig:time}), which 
 explains the increase of the $\mathrm{SNR}_\mathrm{tot}$. A combination of these arguments leads 
 to the peak of the $\mathrm{SNR}_\mathrm{tot}$ at an intermediate $\rho_{\mathrm{p}0}$ value seen
 in Figure \ref{fig:SNRtot_rho}. However, note that in addition to neglecting the 
 initial hyperbolic encounter and the final coalescence/ringdown segments of the signal, our waveform
 model also neglects contributions of spherical moments beyond the quadrupole-order and deviations from
 a precessing Keplerian orbit \citep{Davisetal1972,Bertietal2010,Healyetal2016}. This approximation may
 not be valid for low $\rho_{\mathrm{p}0}$ (particularly for $\rho_{\mathrm{p}0} \leq 10 $). The $\mathrm{SNR}
 _\mathrm{tot}$ is expected to be underestimated in this region in Figure \ref{fig:Delta_rho}.
 

\begin{table*}
\centering  
   \begin{tabular}{@{}cccccccccc}
     \hline
      $\rho_\mathrm{p0}$ & $10$ &$10$& $10$ & $20$ & $20$ & $20$ & Circular& Circular& Circular \\
      quantile & $10\%$ & $50\%$ & $90\%$ & $10\%$ & $50\%$ & $90\%$ &$10\%$ & $50\%$ & $90\%$ \\
     \hline\hline
   $\Delta ( \mathrm{ln} D_\mathrm{L})$  &  $1.13\,(-2)$  &  $2.49\,(-2)$  &  $7.21\,(-2)$  &  $1.89\,(-2)$  &  $4.36\,(-2)$  &  $0.154$  &  $8.81\,(-3)$  &  $5.86\,(-2)$  &  $0.311$  \\
   $\Delta\Omega_N \,\, [\mathrm{sr}]$  &  $1.13\,(-5)$  &  $1.46\,(-4)$  &  $7.23\,(-4)$  &  $2.21\,(-4)$  &  $7.28\,(-4)$  &  $2.89\,(-3)$  &  $5.65\,(-5)$  &  $1.58\,(-3)$  &  $6.67\,(-3)$  \\
   $\Delta\Omega_L \,\, [\mathrm{sr}]$  &  $7.89\,(-5)$  &  $3.85\,(-3)$  &  $0.13$  &  $2.21\,(-5)$  &  $1.13\,(-2)$  &  $0.85$  &  $4.36\,(-4)$  &  $2.21\,(-2)$  &  $9.6$  \\
   $a_N \,\, [\mathrm{deg}]$  &  $0.14$  &  $0.62$  &  $2.42$  &  $0.57$  &  $1.36$  &  $4.07$  &  $0.28$  &  $2.01$  &  $5.16$  \\
   $b_N \,\, [\mathrm{deg}]$  &  $0.11$  &  $0.22$  &  $0.48$  &  $0.26$  &  $0.52$  &  $0.98$  &  $0.41$  &  $0.76$  &  $1.44$  \\
   $a_L \,\, [\mathrm{deg}]$  &  $0.96$  &  $2.78$  &  $20.04$  &  $1.44$  &  $4.65$  &  $45.48$  &  $0.78$  &  $6.98$  &  $2.02\,(+2)$  \\
   $b_L \,\, [\mathrm{deg}]$  &  $0.55$  &  $1.33$  &  $7.51$  &  $0.84$  &  $2.46$  &  $20.58$  &  $0.32$  &  $3.41$  &  $50.63$  \\
   $\Delta \Phi_c \,\, [\mathrm{rad}]$  &  $9.27\,(-2)$ & $0.23$ & $0.71$  &  $0.29$  &  $0.66$  &  $2.44$  &  $0.36$  &  $0.72$  &  $57.91$  \\
   $\Delta t_c \,\, [\mathrm{ms}]$  &  $4.32\,(-2)$ & $8.41\,(-2)$ & $0.181$  &  $9.28\,(-2)$  & $0.167$  & $0.311$ & $9.41\,(-2)$ & $1.582$ & $3.011$  \\
   $\Delta ( \mathrm{ln} \mathcal{M}_z )$  &  $3.53\,(-5)$  &  $6.17\,(-5)$  &  $1.71\,(-4)$  &  $1.04\,(-5)$  &  $1.81\,(-5)$  &  $3.42\,(-5)$  &  $1.36\,(-3)$  &  $2.34\,(-3)$  &  $4.27\,(-3)$  \\
   $\Delta (\mathrm{ln} M_{{\rm tot},z})$  &  $5.42\,(-4)$  &  $9.51\,(-4)$  &  $2.43\,(-3)$  &  $2.82\,(-4)$  &  $4.81\,(-4)$  &  $9.18\,(-4)$  &  $5.88\,(-3)$  &  $1.13\,(-2)$  &  $1.81\,(-2)$  \\
   $\Delta e_\mathrm{LSO}$  &  $1.18\,(-4)$  &  $2.16\,(-4)$  &  $5.83\,(-4)$  &  $2.39\,(-5)$  &  $3.19\,(-5)$  &  $5.88\,(-5)$  &  $-$  &  $-$  &  $-$  \\
   $\Delta e_0$  &  $1.44\,(-3)$  &  $2.16\,(-3)$  &  $3.95\,(-3)$  &  $1.72\,(-3)$  &  $2.91\,(-3)$  &  $5.79\,(-3)$  &  $-$  &  $-$  &  $-$  \\
   $\Delta \rho_{\mathrm{p}0}$  &  $6.64\,(-3)$  &  $1.08\,(-2)$  &  $2.28\,(-2)$  &  $1.33\,(-2)$  &  $2.29\,(-2)$  &  $4.58\,(-2)$  &  $-$  &  $-$  &  $-$  \\
   $\Delta \gamma_c \,\, [\mathrm{rad}]$  &  $0.04$  &  $0.11$  &  $0.47$  &  $0.14$  &  $0.34$  &  $1.51$  &  $-$  &  $-$  &  $-$  \\
\hline\hline
   \end{tabular} 
   \caption{ \label{tab:PrecBin} The $10\%$, $50\%$, and $90\%$ quantile of measurement errors for parameters
   of $30 \, \mathrm{M}_{\odot}-30 \, \mathrm{M}_{\odot}$ precessing eccentric BH binaries with initial eccentricity $e_0=0.9$ and
   dimensionless pericenter distance $\rho_{\mathrm{p}0}=10$ and $20$, and circular binaries at distance $D_
   \mathrm{L}=100 \, \mathrm{Mpc}$, random source sky location and orientation using the detector network 
   in Table \ref{tab:DetCord}. Here $(\Omega_N, a_N, b_N)$ are respectively the area, semi-major and semi-minor 
   axes of the 2D error ellipse corresponding to the source's sky direction, and similarly for $(\Omega_L, 
   a_L, b_L)$ describing the source's orbital plane orientation (i.e. angular momentum vector direction). 
   Note that $e_{\rm LSO}=0.187$ and $0.059$ for $\rho_{\mathrm{p}0}=10$ and $20$, respectively, if $e_0=0.9$
   is assumed. In the circular limit the binary forms outside of the sensitive frequency band of the detector 
   network, and $\Delta e_0\rightarrow \infty$ and $\Delta \rho_{\mathrm{p}0}\rightarrow \infty$. We adopt the 
   following notation in the table: $1.13\,(-2)=1.13 \times 10^{-2}$. }
\end{table*}

\begin{table*}
\centering  
   \begin{tabular}{@{}ccccccc}
     \hline
      $\rho_{\mathrm{p}0}$ & $10$ & $10$ & $10$ & $20$ & $20$ & $20$  \\
      $e_0$ & $0.9$ &$0.95$& $0.97$ & $0.9$ & $0.95$ & $0.97$  \\
      $e_\mathrm{LSO}$  &  $0.1872$  &  $0.1932$  &  $0.1956$  &  $5.89\,(-2)$   &  $6.07\,(-2)$ &  $6.15\,(-2)$ \\
      $\mathrm{SNR}_\mathrm{tot}$  &  $251.1$  &  $257.7$  &  $260.4$  &  $179.4$   &  $182.1$ &  $182.9$ \\
     \hline\hline
   $\Delta ( \mathrm{ln} D_\mathrm{L})$  &  $3.92\,(-2)$  &  $3.79\,(-2)$ &  $3.76\,(-2)$  &  $6.41\,(-2)$  
   &  $6.31\,(-2)$   &  $6.27\,(-2)$  \\
   $\Delta\Omega_N \,\, [\mathrm{sr}]$  &  $8.37\,(-5)$  &  $7.78\,(-5)$  &  $7.53\,(-5)$  &  $4.06\,(-4)$ 
   &  $3.91\,(-4)$   &  $3.85\,(-4)$  \\
   $\Delta\Omega_L \,\, [\mathrm{sr}]$  &  $9.68\,(-3)$  &  $9.07\,(-3)$  &  $8.84\,(-3)$  &  $2.74\,(-2)$ 
   &  $2.65\,(-2)$  &  $2.61\,(-2)$  \\
   $a_N \,\, [\mathrm{deg}]$  &  $0.501$  &  $0.484$  &  $0.477$  &  $0.994$  &  $0.977$  &  $0.970$  \\
   $b_N \,\, [\mathrm{deg}]$  &  $0.174$  &  $0.168$  &  $0.165$  &  $0.427$  &  $0.418$  &  $0.415$  \\
   $a_L \,\, [\mathrm{deg}]$  &  $4.31$  &  $4.18$  &  $4.13$    &  $6.96$ &  $6.85$ &  $6.81$  \\
   $b_L \,\, [\mathrm{deg}]$  &  $2.35$  &  $2.26$  &  $2.23$  &  $4.12$  &  $4.04$  &  $4.01$  \\   
   $\Delta \Phi_c \,\, [\mathrm{rad}]$  &  $0.331$ & $0.313$  &  $0.308$  & $1.102$ &  $1.049$  &  $1.031$  \\
   $\Delta t_c \,\, [\mathrm{ms}]$  &  $9.68\,(-2)$ & $9.08\,(-2)$  & $8.88\,(-2)$  & $0.198$ &  $0.192$ 
   & $0.189$   \\
   $\Delta ( \mathrm{ln} \mathcal{M}_z )$  &  $4.61\,(-5)$  &  $1.89\,(-5)$  &  $1.39\,(-5)$  &  $1.38\,(-5)$ 
   &  $6.03\,(-6)$  &  $4.33\,(-6)$  \\
   $\Delta (\mathrm{ln} M_{{\rm tot},z})$  &  $7.09\,(-4)$  &  $5.62\,(-4)$  &  $4.85\,(-4)$  &  $3.67\,(-4)$ 
   &  $2.87\,(-4)$  &  $2.38\,(-4)$  \\
   $\Delta e_\mathrm{LSO}$  &  $1.61\,(-4)$  &  $1.28\,(-4)$ &  $1.11\,(-4)$ &  $2.43\,(-5)$  &  $1.93\,(-5)$ 
   &  $1.62\,(-5)$   \\
   $\Delta e_0$  &  $1.95\,(-3)$  &  $1.74\,(-3)$  &  $1.69\,(-3)$ &  $2.25\,(-3)$  &  $2.17\,(-3)$  &  $2.14\,(-3)$   \\
   $\Delta \rho_{\rm p0}$  &  $8.33\,(-3)$  &  $7.45\,(-3)$  & $7.11\,(-3)$  &  $1.71\,(-2)$  &  $1.61\,(-2)$ 
   &  $1.56\,(-2)$   \\
   $\Delta \gamma_c \,\, [\mathrm{rad}]$  &  $0.169$  &  $0.162$ &  $0.159$  &  $0.556$   &  $0.529$  &  $0.520$    \\
\hline\hline
   \end{tabular}  
   \caption{ \label{tab:PrecBine0rho0Dep} Measurement errors for parameters of $30 \, \mathrm{M}_{\odot}-30 \, \mathrm{M}_{\odot}$ 
   precessing eccentric BH binaries with initial eccentricities $e_0=0.9$, $0.95$, and $0.97$ for initial
   dimensionless pericenter $\rho_{p0}=10$ and $20$, luminosity distance $D_\mathrm{L}=100 \, \mathrm{Mpc}$, 
   and arbitrarily fixed source direction $(\theta_N,\phi_N)=(\pi/2,\pi/3)$ and orientation $(\theta_L,
   \phi_L)=(\pi/4,\pi/5)$ for the detector network specified in Table \ref{tab:DetCord}. }
\end{table*} 
 
 
 
\subsection{Parameter measurement errors} 
\label{subsec:ResParamDist} 
 
\begin{figure*}
    \centering 
    \includegraphics[width=79mm]{Dist_P_lnD} 
    \includegraphics[width=79mm]{Dist_P_lnM}    
    \\
    \includegraphics[width=79mm]{Dist_P_aN}   
    \includegraphics[width=79mm]{Dist_P_bN}
    \\
    \includegraphics[width=79mm]{Dist_P_aL}   
    \includegraphics[width=79mm]{Dist_P_bL}
\caption{ Smoothed distribution of the measurement errors of parameters measured for $30 \, \mathrm{M}_{\odot}-30 \,
 \mathrm{M}_{\odot}$ precessing eccentric BH binaries with initial eccentricity $e_0 = 0.9$ and dimensionless pericenter 
 distance $\rho_{\mathrm{p}0}=10$ and $20$ (green dashed and red dash-dot) and for similar binaries in 
 the circular limit (blue solid) at luminosity distance $D_\mathrm{L} = 100 \, \mathrm{Mpc}$ with a
 random source direction and binary orientation. We have assumed the detector network specified in 
 Table \ref{tab:DetCord}. \emph{Top row:} Distribution of the relative measurement error of luminosity 
 distance, $\Delta D_{\rm L} / D_{\rm L} = \Delta (\ln D_{\rm L})$, and redshifted chirp mass, $\Delta 
 \mathcal{M}_{z} / \mathcal{M}_{z} = \Delta (\ln \mathcal{M}_{z})$. \emph{Middle row:} Distribution of 
 semi-major axis of the sky position error ellipse, $a_N$, and its semi-minor axis, $b_N$. \emph{Bottom 
 row:} Distribution of the semi-major axis of the binary's orbital plane orientation error ellipse, $a_L$,
 and its semi-minor axis, $b_L$.  } \label{fig:Dist1}
\end{figure*} 


\begin{figure*}
    \centering 
    \includegraphics[width=79mm]{Dist_P_e0}   
    \includegraphics[width=79mm]{Dist_P_rhop0}
    \\
    \includegraphics[width=79mm]{Dist_P_eLSO}
\caption{ Smoothed distribution of the measurement errors of various parameters measured only for 
 precessing eccentric binaries. Similar to Figure \ref{fig:Dist1}, these distributions correspond
 to a Monte Carlo sample of $4500$ binaries with random source direction and binary orientation for 
 $30 \, \mathrm{M}_{\odot} - 30 \, \mathrm{M}_{\odot}$ precessing eccentric BH binaries with initial eccentricity $e_0 = 
 0.9$ and dimensionless pericenter distance $\rho_{\mathrm{p}0} = 10$ and $20$ (green dashed and red 
 dash-dot). The source distance is also fixed at $D_\mathrm{L} = 100\, \mathrm{Mpc}$, and the detector 
 network is specified in Table \ref{tab:DetCord}. \emph{Top left:} Distribution of the measurement 
 error in initial orbital eccentricity, $e_0$. \emph{Top right:} Distribution of the measurement error
 in initial dimensionless pericenter distance, $\rho_{\mathrm{p}0}$. \emph{Bottom:} Distribution of the 
 measurement error in the eccentricity at the last stable orbits, $e_{\rm LSO}$. }  \label{fig:Dist2}
\end{figure*} 

 We present the measurement accuracy for the final eccentricity at the LSO for parameters grouped as  
\begin{align} 
  \bm{\lambda}_\mathrm{slow} & = \lbrace \mathrm{ln} (D_\mathrm{L}), \theta_N, \phi_N, \theta_L, 
  \phi_L \rbrace \, ,  \\
  \bm{\lambda}_\mathrm{fast} & = \lbrace \Phi_c,  t_c, \mathrm{ln} (\mathcal{M}_z), 
  \mathrm{ln} (M_{{\rm tot},z}), e_\mathrm{LSO}, \gamma_c \rbrace \,\
\end{align} 
 \citep{Kocsisetal2007}. The $\bm{\lambda}_\mathrm{fast}$ fast parameters are related to the high 
 frequency GW phase, while the $\bm{\lambda}_\mathrm{slow}$ slow parameters appear only in the 
 slowly-varying amplitude of the GW signal. Slow parameters are mostly determined by a comparison 
 of the GW signals measured by the different detectors in the network. For the polar angles $(\theta_N,
 \phi_N)$ describing the source direction, we calculate the minor and major axes ($a_N,b_N$) of the 
 corresponding 2D sky location error ellipse and its area ($\Omega_N= \pi a_N b_N$), and we do the 
 same for the binary orientation error ellipse $(a_L,b_L)$ and its area ($\Omega_N=\pi a_L b_L$).

 Figures \ref{fig:Dist1} and \ref{fig:Dist2} show the distribution of the measurement errors for randomly
 chosen source sky position and binary orientation for $\rho_{\mathrm{p}0}=10$ and $20$, and for the circular
 limit $\rho_{\mathrm{p}0} \rightarrow \infty$ (see Appendix \ref{sec:circlim}), while Table \ref{tab:PrecBin}
 shows the $10 \%$, $50\%$, and $90\%$ quantiles of the error distributions. Compared to a $\rho_{\mathrm{p}0}
 = 20$ binary, a $\rho_{\mathrm{p}0} = 10$ binary is more eccentric throughout its evolution, which leads to
 a higher $\mathrm{SNR}_\mathrm{tot}$, and most of its measurement errors are smaller. There are, however, 
 exceptions to this finding: the fast parameters such as the mass parameters and the eccentricity have higher errors for $\rho_{\mathrm{p}0}
 = 10$ than for $\rho_{\mathrm{p}0} = 20$ (see discussion below). 
  
 Many of the binaries in galactic nuclei\footnote{particularly the heavy BHs therein \citep{Gondanetal2017}} form with very high $e_0$, close to unity, in single-single
 encounters due to GW emissions \citep{OLearyetal2009,Gondanetal2017}. However, similarly to the finding that 
 the $\mathrm{SNR}_\mathrm{tot}$ does not increase significantly for $1 > e_0 \geq 0.9$, we find that $\bm{\lambda}
 _\mathrm{slow}$ parameter errors do not improve due to the early very eccentric evolutionary period beyond $e > 
 0.9$ (repeated burst phase) compared to waveforms with $e_0=0.9$ as shown in Table \ref{tab:PrecBine0rho0Dep}. 
 However, some $\bm{\lambda}_\mathrm{fast}$ parameters' measurement errors improve more 
 significantly with $e_0 > 0.9$. In particular, the measurement errors of the mass parameters ($\mathcal{M}_z,
 M_{{\rm tot},z}$) improve by a factor of $\sim 2$, and the measurement error of $e_\mathrm{LSO}$ improves 
 by $\sim 50\%$ if increasing $e_0$ from $0.9$ for $0.97$. This difference is due to the fact that eccentricity
 modifies the GW phase significantly, which affects the determination of $\bm{\lambda}_\mathrm{fast}$ parameters
 only.

\begin{figure*}
    \centering 
    \includegraphics[width=79mm]{DeltaDL_rho}  
    \includegraphics[width=79mm]{DeltaMchp_rho} \\
    \includegraphics[width=79mm]{aNbN_rho}  
    \includegraphics[width=79mm]{aLbL_rho} \\
    \includegraphics[width=79mm]{e0_rho}  
    \includegraphics[width=79mm]{rhop0_rho} \\
    \includegraphics[width=79mm]{eLSO_rho}
\caption{ The measurement error of source parameters as a function of initial dimensionless pericenter 
 distance with all other binary parameters fixed as in Figure \ref{fig:SNRtot_rho}. \emph{First row left:}
 luminosity distance $\Delta D_\mathrm{L}/ D_\mathrm{L} = \Delta (\mathrm{ln} D_\mathrm{L})$. 
 \emph{First row right:} redshifted chirp mass, $\Delta \mathcal{M}_z/\mathcal{M}_z= \Delta 
 (\mathrm{ln} \mathcal{M}_z)$. \emph{Second row left:} Semi-major and semi-minor axes of the sky 
 position error ellipse $a_N$ and $b_N$. \emph{Second row right:} Semi-major and semi-minor axes of
 the error ellipse for the binary orbital plane normal vector direction, $a_L$ and $b_L$. \emph{Third 
 row left:} Initial orbital eccentricity, $\Delta e_0$. \emph{Third row right:} Initial dimensionless 
 pericenter distance, $\Delta \rho_{\rm p0}$. \emph{Fourth row:} Eccentricity at the last stable orbit, 
 $\Delta e_{\rm LSO}$. The measurement error of parameters ($\mathrm{ln} (D_\mathrm{L}), \mathrm{ln} 
 (\mathcal{M}_z), a_N, b_N, a_L, b_L$) converge asymptotically to the value of precessing eccentric
 binaries in the circular limit for high $\rho_\mathrm{p0}$, and the measurement error of parameters 
 ($\Delta e_0, \Delta \rho_{\mathrm{p}0}, \Delta e_{\rm LSO}$) increase rapidly with $\rho_{\mathrm{p}0}$
 for high $\rho_{\mathrm{p}0}$. We find similar trends with $\rho_{\mathrm{p}0}$ for other random 
 choices of binary direction and orientation (not shown). Note that the measurement error of $\Delta 
 e_{\rm LSO}$ is undetermined for high $\rho_{\mathrm{p}0}$ because the Fisher matrix algorithm becomes
 invalid for this parameter in this regime, see Section \ref{subsec:ResParamDist} for details. }
 \label{fig:Delta_rho}
\end{figure*} 

 We calculate the parameter measurement errors for precessing highly eccentric BH binaries as a 
 function of $\rho_{\mathrm{p}0}$ for some arbitrarily fixed binary direction and orientations. 
 For one such binary direction and orientation Figure \ref{fig:Delta_rho} shows the 
 $\rho_{\mathrm{p}0}$ dependence of $\Delta D_{\rm L} / D_{\rm L}$, $\Delta \mathcal{M}_z / 
 \mathcal{M}_z$, $\Delta e_0$, $\Delta \rho_ {\mathrm{p}0}$, $\Delta e_{\rm LSO}$, semi-major 
 and semi-minor axes of the sky position error ellipse $(a_N,b_N)$, and semi-major and semi-minor
 axes of the error ellipse for the binary orbital plane normal vector direction $(a_L,b_L)$. Not
 e that we find similar trends with $\rho_{\mathrm{p}0}$ for other random choices of binary direction
 and orientation. We find that measurement errors systematically decrease with decreasing $\rho_
 {\mathrm{p}0}$ for precessing highly eccentric binaries relative to similar binaries in the circular
 limit (Appendix \ref{sec:circlim}), and the errors have a minimum in the range $8 < \rho_{\mathrm{p}0} 
 < 80$ and deteriorate rapidly for $\rho_{\mathrm{p}0} < 8$. The latter is due to the rapid decrease of 
 the $\mathrm{SNR}_\mathrm{tot}$ in that range. The $\rho_{\mathrm{p}0}$ dependence of 
 $\Delta D_\mathrm{L} / D_\mathrm{L}$ and the principal axes of the sky position and binary orientation
 error ellipses $(a_N,b_N)$ and $(a_L,b_L)$ (i.e. quantities derived from slow parameters, see Section 
 \ref{sec:ParamSpace}) are qualitatively similar to that of $1/\mathrm{SNR}_\mathrm{tot}$ in the complete
 range of $\rho_{\mathrm{p}0}$ (i.e. they decrease rapidly with $\rho_{\mathrm{p}0}$ for low $\rho_
 {\mathrm{p}0}$, have a minimum at moderate $\rho_{\mathrm{p}0}$, and converge asymptotically to the 
 value of precessing highly eccentric binaries in the circular limit), see Figure \ref{fig:ecc_per_circ} 
 for details. However, Figure \ref{fig:Delta_rho} shows that the chirp mass errors have a minimum at 
 much higher $\rho_{\mathrm{p}0}$, i.e. between $50 - 60$ and $20 - 40$ for $10 \,\mathrm{M}_{\odot} - 10 \, \mathrm{M}_{\odot}$ 
 and $30 \, \mathrm{M}_{\odot} -30 \, \mathrm{M}_{\odot}$ precessing highly eccentric BH binaries, respectively. The main reason 
 for the different behavior of the chirp mass from the distance and angular errors is the fact that the 
 chirp mass is a fast parameter, while the distance and angular parameters are slow parameters. Slow 
 parameters are insensitive to the GW phase perturbations, and depend on the GW amplitude, which is set
 by the $\mathrm{SNR}_\mathrm{tot}$. The SNR of the early part of the waveform near the low-frequency 
 noise wall of the detector is small. However, fast parameters depend sensitively on the GW phase, and the 
 GW phase accumulates mostly at low frequencies, since the residence time (i.e. $\nu/\dot \nu$) is largest
 at low orbital frequencies. Thus, the fast parameters' errors are minimized for binaries which form with 
 $e_0 \sim 1$ with a $\rho_{\rm p0}$ value for which the GW characteristic frequency is near the detectors' 
 minimum frequency. The peak of the spectrum is initially at $f_{\min}$ if $\rho_{\rm p0} = 40.9 [M_{\rm tot}
 /(20 \, \mathrm{M}_{\odot})]^{-2/3} (f_{\min} / 10\,{\rm Hz})^{-2/3}$ \citep{Gondanetal2017}. A slightly lower value of 
 $f_{\min} \sim 7 \, {\rm Hz}$ leads to values that represent the minimum of the fast parameters. This also 
 leads to the result observed in Figure \ref{fig:Dist1} that these parameters have higher errors for $\rho_
 {\mathrm{p}0} = 10$ than for $\rho_{\mathrm{p}0} = 20$. 
 
\begin{figure}
    \centering 
    \includegraphics[width=78mm]{ecc_per_circ_MA30_MB30}  
\caption{ The initial dimensionless pericenter distance dependence of $\Delta D_\mathrm{L} / 
 D_\mathrm{L}$, the principal axes of the sky localization and binary orientation error ellipses 
 $(a_N,b_N)$ and $(a_L,b_L)$, and $1 / \mathrm{SNR}_\mathrm{tot}$ for $30 \, \mathrm{M}_{\odot}-30 \, \mathrm{M}_{\odot}$ 
 precessing highly eccentric BH binaries with all other binary parameters fixed as in Figure 
 \ref{fig:SNRtot_rho}. The values are with respect to those of circular binaries and match those 
 in Figure \ref{fig:Delta_rho}. Note that we find similar trends with $\rho_{\mathrm{p}0}$ and 
 $e_0$ for other random choices of binary direction and orientation (not shown). }
 \label{fig:ecc_per_circ}
\end{figure}  
 
 In Figure \ref{fig:Delta_rho}, note that $\Delta e_{\mathrm{LSO}}$ errors are relatively small for 
 relatively high $\rho_{\mathrm{p}0}$ up to $\rho_{\mathrm{p}0} \sim 200$. At high $\rho_{\mathrm{p}0}$, 
 the orbital eccentricity approaches zero when it enters the aLIGO band, and $\Delta e_{\mathrm{LSO}}$ 
 increases. We note that the posterior probability distribution function of $e_{\mathrm{LSO}}$ is 
 well-defined even in the circular limit $\rho_{\mathrm{p}0} \rightarrow \infty$, and $\Delta e_
 {\mathrm{LSO}}$ is finite for a given confidence region. However, the Fisher matrix algorithm becomes
 invalid in this regime as the signal is not approximated well by its linear Taylor expansion with respect
 to the $\Delta e_{\mathrm{LSO}}$ parameter, since its first $e_{\mathrm{LSO}}$ derivative vanishes in the 
 circular limit. Therefore the true asymptotic value of $\Delta e_{\mathrm{LSO}}$ for high $\rho_{\mathrm{p}0}$
 cannot be recovered with the Fisher matrix technique used in this paper. Further, note that $\Delta e_0$ 
 and $\Delta \rho_{\mathrm{p}0}$ also increase rapidly with $\rho_{\mathrm{p}0}$ for high $\rho_{\mathrm{p}0}$.
 This is due to the fact that for these parameters the binary forms with a pericenter frequency smaller than 
 the minimum frequency of the detector network, and the information on $e_0$ and $\rho_{\mathrm{p}0}$ is 
 limited to higher harmonics with small power. Thus, these parameters indeed have a very high error and 
 become indeterminate in the circular limit. The fact that the relative error of $e_0$ and $\rho_{\mathrm{p}0}$ 
 can be less than $\sim 5\%$ percent in the range $5 < \rho_{\mathrm{p}0} < 50$ ($ 6 < \rho_{\mathrm{p}0} < 
 100$) for $30 \, \mathrm{M}_{\odot} - 30 \, \mathrm{M}_{\odot}$ ($10 \, \mathrm{M}_{\odot} - 10 \, \mathrm{M}_{\odot}$) precessing highly eccentric BH binaries 
 implies that the GW detections might have the potential to constrain the formation environment of these system 
 \citep{OLearyetal2009,Cholisetal2016,Rodriguezetal2016a,Rodriguezetal2016b,Chatterjeeetal2017,Gondanetal2017,Kocsisetal2017,SamsingRamirezRuiz2017,SilsbeeTremaine2017}. 
 
 Furthermore, we found from numerical investigations that $\Delta e_0$ does not correlate
 significantly with other parameters' errors, which is due to the fact that $e_0$ is measured from 
 the truncation of the signal for $e>e_0$ at the start of the waveform, while other parameters of 
 a precessing eccentric binary are measured from the inspiral rate (Section \ref{sec:ParamSpace}). 
 However, $\Delta \rho_{\mathrm{p}0}$ behaves differently from $\Delta e_0$ in this regard, which 
 is due to the fact that $\rho_{\mathrm{p}0}$ is determined by $e_{\mathrm{LSO}}$ in Equation 
 (\ref{eq:deltarhop0}), and $e_{\mathrm{LSO}}$ depends on the mass parameters. 




\subsection{Comparison with previous results}
\label{sebsec:CompWithPrevRes}
 
 In this paper, we have determined the SNR and the expected accuracy with which the 
 aLIGO-AdV-KAGRA detector network may determine the parameters that describe highly eccentric BH 
 binaries, and investigated how these quantities depend on the initial pericenter distance $\rho_
 {\mathrm{p}0}$ and initial eccentricity $e_0$. There are some previous studies that also made 
 similar investigations for eccentric compact binaries with significant differences 
 \citep{Yunesetal2009,KyutokuSeto2014,Sunetal2015,Maetal2017}. They considered different detector 
 networks, applied different waveform models, and used different definitions for $e_0$ and $\rho_
 {\mathrm{p}0}$. As a consequence, only a qualitative comparison is possible with those results, 
 which we discuss in this section. At the end of this section, we compare our results for the 
 measurement errors in the circular limit with those presented in previous studies.
 
 We first compare our results with a previous study for the $\rho_{\mathrm{p}0}$ dependence 
 of the $\mathrm{SNR}_\mathrm{tot}$. Our result for the $\rho_{\mathrm{p}0}$ dependence of the $\mathrm{SNR}_
 \mathrm{tot}$ (Figure \ref{fig:SNRtot_rho}) is qualitatively in agreement with the result of Figure
 2 in \citet{KyutokuSeto2014}, i.e. the $\mathrm{SNR}_\mathrm{tot}$ increases rapidly with $\rho_{\mathrm{p}0}$
 for low $\rho_{\mathrm{p}0}$, peaks at a moderate $\rho_{\mathrm{p}0}$, and converges asymptotically 
 to the value of highly eccentric binaries in the circular limit for high $\rho_{\mathrm{p}0}$. 
       
 In order to compare our results for the $e_0$ dependence of the $\mathrm{SNR}_\mathrm{tot}$ with
 \citet{Yunesetal2009} and \citet{Sunetal2015}, we also set the lower bound of advanced GW detectors' 
 sensitive frequency band to $f_{\min} =20 \, \mathrm{Hz}$. We define $e_{20 \, \mathrm{Hz}}$ to be the eccentricity
 at which the peak GW frequency of the binary defined in \citet{Wen2003} is $f_\mathrm{GW} = 20 \, 
 \mathrm{Hz}$, and evaluate $\rho_{\mathrm{p}0}$ corresponding to $e_0=e_{\mathrm{20\mathrm{Hz}}}$ from Equation
 (37) in \citet{Wen2003} as
 \begin{equation}  \label{eq:rho_GW}
   \rho_{20\mathrm{Hz}} = \left[ (1+e_{20\mathrm{Hz}})^{0.3046} f_\mathrm{GW} \pi M_{{\rm tot},z} \right]^{-2/3} \, .
 \end{equation} 
 We recalculate the distribution of the $\mathrm{SNR}_\mathrm{tot}$ for $10 \, \mathrm{M}_{\odot}-10 \, \mathrm{M}_{\odot}$ 
 precessing eccentric compact binaries with $e_{20\mathrm{Hz}}= (0.1, 0.2, 0.3, 0.4)$. The top panel of Figure
 \ref{fig:SNRtot_Comps} shows that the $\mathrm{SNR}_\mathrm{tot}$ is roughly the same for different $e_
 {20\mathrm{Hz}}$, which is consistent with results presented in Figure 2 in \citet{Sunetal2015} and in the 
 left panel of Figure 8 in \citet{Yunesetal2009}. Moreover, we find that the $\mathrm{SNR}_\mathrm{tot}$
 increases weakly with $e_\mathrm{20\mathrm{Hz}}$, which is in agreement with results in the left panel of 
 Figure 8 in \citet{Yunesetal2009}. Note that this result disagrees with Table 5 in \citet{Sunetal2015}. 
 The $\mathrm{SNR}_\mathrm{tot}$ does not depend significantly on $e_{20\mathrm{Hz}}$ in the range of $[0.1,0.4]$ 
 for $10 \, \mathrm{M}_{\odot}-10 \, \mathrm{M}_{\odot}$ for $\rho_{p0}=\rho_{20\mathrm{Hz}}\sim 28$ as seen in the bottom panel of Figure
 \ref{fig:SNRtot_Comps}.\footnote{For binaries with relatively high $\rho_{\mathrm{p}0}$, binaries are 
 well-circularized by the time their peak GW frequency enters the sensitive frequency band of advanced
 ground-based GW detectors, thus the information about the initial eccentricity vanishes from the detectable
 part of the waveform. This explains the very weak $e_0$ dependence of the $\mathrm{SNR}_\mathrm{tot}$ 
 and of parameter measurement errors for high $\rho_{\mathrm{p}0}$ in the bottom panel of Figure 
 \ref{fig:SNRtot_Comps} and in Figure \ref{fig:Err_LowEccDep}. } For $10 \, \mathrm{M}_{\odot} - 10 \, \mathrm{M}_{\odot}$ binaries
 with $e_0<0.4$ and $\rho_{\mathrm{p}0} \sim 28$, we find that $\rho_{20\mathrm{Hz}}$ is high enough to fall into
 the range of $\rho_{\mathrm{p}0}$ where the $\mathrm{SNR}_\mathrm{tot}$ depends on $e_0$ at the $\sim 10\%$
 level for $e_{20\mathrm{Hz}}<0.4$. The influence of $e_{20\mathrm{Hz}}$ on the $\mathrm{SNR}_\mathrm{tot}$ increases with
 $M_\mathrm{tot}$ since in this case $\rho_{20\mathrm{Hz}}$ is lower as shown by Equation (\ref{eq:rho_GW}). Thus, 
 the influence of $e_{20\mathrm{Hz}}$ on the distribution of $\mathrm{SNR}_\mathrm{tot}$ is more significant for 
 higher-mass low-eccentricity binaries. Examples for this characteristic of the $\mathrm{SNR}_\mathrm{tot}$ 
 are seen in Figure 8 in \citet{Yunesetal2009}.


\begin{figure}
    \center
    \includegraphics[width=78mm]{SNRtot_Comp}  \\
    \includegraphics[width=78mm]{SNRtot_e0dep}
\caption{ 
     \emph{Top panel}: Smoothed probability density function of the $\mathrm{SNR}_\mathrm{tot}$
     for $10 \, \mathrm{M}_{\odot}-10 \, \mathrm{M}_{\odot}$ precessing eccentric BH binaries. In these calculations we set the 
     lower bound of the detectors' sensitive frequency band to $20 \, \mathrm{Hz}$. Here $e_{20\mathrm{Hz}}$ 
     and $\rho_{20\mathrm{Hz}}$ represent the initial orbital eccentricity and initial dimensionless pericenter 
     distance at which the peak GW frequency \citep{Wen2003} of the binary is $f_\mathrm{GW} = 20 \, 
     \mathrm{Hz}$. Other details of the calculations are the same as in Figure \ref{fig:SNR_Prec}. 
     \emph{Bottom panel:} The same as in Figure \ref{fig:SNRtot_rho} but for $10 \, \mathrm{M}_{\odot} -10 \, \mathrm{M}_{\odot}$
     precessing eccentric BH binaries with initial eccentricities $e_0= (0.1, 0.2, 0.3, 0.4, 0.9)$ as a 
     function of $\rho_{\mathrm{p}0}$. We also represent the $e_0 = 0.9$ curve in this plot in order to 
     show that $\mathrm{SNR}_\mathrm{tot}$ is a strictly monotonically increasing function of $e_0$ over 
     the full range of $\rho_{\mathrm{p}0}$. We find similar trends with $\rho_{\mathrm{p}0}$ 
     and $e_0$ for other random choices of binary direction and orientation (not shown). }   
     \label{fig:SNRtot_Comps}
\end{figure} 


\begin{figure*}
    \centering 
    \includegraphics[width=78mm]{lnDL_Err_e0Dep} 
    \includegraphics[width=78mm]{lnMchpz_Err_e0Dep} \\
    \includegraphics[width=78mm]{aN_deg_ErrDep} 
    \includegraphics[width=78mm]{aL_deg_ErrDep} 
\caption{ The initial pericenter distance and eccentricity dependence of the measurement errors
 of various parameters for $10 \, \mathrm{M}_{\odot}-10 \, \mathrm{M}_{\odot}$ precessing eccentric BH binaries for a specific 
 sky position and inclination chosen as in Figure \ref{fig:SNRtot_rho}. Different lines show different 
 initial eccentricities as labeled. \emph{Top left:} luminosity distance $\Delta D_\mathrm{L}/ D_\mathrm{L} 
 = \Delta (\mathrm{ln} D_\mathrm{L})$. \emph{Top right:} redshifted chirp mass, $\Delta \mathcal{M}_z/
 \mathcal{M}_z= \Delta (\mathrm{ln} \mathcal{M}_z)$. \emph{Bottom left:} Semi-major axis of the sky 
 position error ellipse, $a_N$. \emph{Bottom right:} Semi-major axis of the error ellipse for the 
 binary orbital plane normal vector direction, $a_L$. We find similar trends with $\rho_{\mathrm{p}0}$ 
 and $e_0$ for other random choices of binary direction and orientation (not shown). We also show
 the $e_0 = 0.9$ curve in this plot to illustrate that measurement errors of binary parameters are a
 strictly monotonically increasing functions of $e_0$ over the full range of $\rho_{\mathrm{p}0}$. }
 \label{fig:Err_LowEccDep}
\end{figure*} 
 
 Finally, we compare our results with previous studies for the $e_0$ dependence of measurement
 errors of parameters describing eccentric binaries. \citealt{Sunetal2015,Maetal2017} set the initial 
 orbital parameters to be $e_{20\mathrm{Hz}}$ ($e_{10\mathrm{Hz}}$) and $\rho_{20\mathrm{Hz}}$ ($\rho_{10\mathrm{Hz}}$). For various 
 values of $e_{20\mathrm{Hz}}$ ($e_{10\mathrm{Hz}}$) in the range $[0.1, 0.2, 0.3, 0.4]$ and the corresponding values
 of $\rho_{20\mathrm{Hz}}$ ($\rho_{10\mathrm{Hz}}$), they determined the measurement accuracies for various parameters
 of eccentric binaries. Since they applied different waveform models and different parameters 
 describing the eccentric binaries, we resort to a qualitative comparison. We repeated the analysis of 
 Figure \ref{fig:Delta_rho} for $e_0= (0.1, 0.2, 0.3, 0.4)$ and determined the measurement error of 
 parameters as a function of $\rho_{\mathrm{p}0}$ as shown in Figure \ref{fig:Err_LowEccDep}. We find 
 qualitative agreement with \citet{Sunetal2015,Maetal2017}. The measurement accuracies of parameters 
 increase strictly monotonically with $e_0$. 
 
 Previous papers have investigated the $M_\mathrm{tot}$ dependency of the measurement errors for $\lbrace 
 t_c, \Phi_c, \mathrm{ln}(\mathcal{M}_{z}), \mathrm{ln}(\eta) \rbrace$ by using different PN order waveform 
 models for non-spinning inspiraling binaries for a fixed SNR in a single aLIGO type detector. Previous 
 results showed that the measurement accuracy of these parameters decreases with increasing $M_\mathrm{tot}$  
 for $2.8 \, \mathrm{M}_{\odot} \leq M_\mathrm{tot} \leq 20 \, \mathrm{M}_{\odot}$ provided that the SNR accumulated in one GW detector 
 is fixed, see Table 1 in \citet{Arunetal2005} and references therein. Therefore, we determined the measurement 
 errors in the circular limit for $\lbrace t_c, \Phi_c, \mathrm{ln}(\mathcal{M}_{z}), \mathrm{ln}(\eta) \rbrace$
 for a qualitative comparison.\footnote{A quantitative agreement is not expected since our precessing waveform 
 model differs from the waveform models in those studies.} To calculate the measurement error of the $\mathrm{ln}
 (\eta)$ parameter we use the fact that $\eta= (\mathcal{M}_z M_{{\rm tot},z}^{-1})^{5/3}$ and so 
\begin{equation}  \label{eq:Delta_eta} 
  \frac{\langle\Delta \eta^2\rangle}{\eta^2}  = \frac{25}{9} \frac{\langle\Delta \mathcal{M}_z^2\rangle}{
  \mathcal{M}_z^2}  +\frac{25}{9} \frac{ \langle\Delta M_{{\rm tot},z}^2\rangle }{ M_\mathrm{tot,z}^2 } 
  + \frac{50}{9} \frac{ \langle \Delta \mathcal{M}_z \Delta M_{{\rm tot},z} \rangle }{ 
  \mathcal{M}_z M_{{\rm tot},z} } \, .
\end{equation}
 In agreement with the 1PN order case in \citet{Arunetal2005}, we find that $\Delta t_c$, $\Delta \Phi_c$,
 $\Delta \mathcal{M}_{z} / \mathcal{M}_z$, and $\Delta \eta / \eta$ increase with $M_\mathrm{tot}$ for fixed 
 $\mathrm{SNR}_\mathrm{tot}$ (Table \ref{Table:CircLimParams}). Such a qualitative agreement is expected since 
 the adopted precessing eccentric waveform approximates the full 1PN waveform in its most important features, and 
 the $M_\mathrm{tot}$-dependent trends of error distributions do not depend on the number of detectors or on the 
 sky position or angular momentum unit vectors of the source. 
 
\begin{table} 
\centering 
   \begin{tabular}{@{}ccccc} 
      $m_A-m_B$  &  $\Delta t_c$  &  $\Delta \Phi_c $  &  $\Delta \mathcal{M}_{z} / \mathcal{M}_z$ 
      &  $\Delta \eta / \eta$  \\ 
     \hline\hline 
   $10 \, \mathrm{M}_{\odot}-10 \, \mathrm{M}_{\odot}$  &  $0.30$  &  $0.32$  &  $3.6 \times 10^{-4}$  &  $3.4 \times 10^{-3}$  \\
   \hline
   $15 \, \mathrm{M}_{\odot}-15 \, \mathrm{M}_{\odot}$  &  $0.55$  &  $0.39$  &  $8.7 \times 10^{-4}$  &  $6.5 \times 10^{-3}$  \\
   \hline
   $20 \, \mathrm{M}_{\odot}-20 \, \mathrm{M}_{\odot}$  &  $0.94$  &  $0.49$  &  $1.5 \times 10^{-3}$  &  $10^{-2}$  \\
   \hline
   $25 \, \mathrm{M}_{\odot}-25 \, \mathrm{M}_{\odot}$  &  $1.38$  &  $0.57$  &  $2.2 \times 10^{-3}$  &  $1.3 \times 10^{-2}$  \\
   \hline
   $30 \, \mathrm{M}_{\odot}-30 \, \mathrm{M}_{\odot}$  &  $1.84$  &  $0.62$  &  $2.8 \times 10^{-3}$  &  $1.5 \times 10^{-2}$  \\
   \hline\hline 
   \end{tabular} 
   \caption{ Errors in $t_c$ (msec), $\Phi_c$ (rad) and for the relative errors
   in $\mathcal{M}_z$ and $ \eta $ in the circular limit for equal-mass binaries 
   for a specific sky position and inclination
   $\theta_N=\pi/2$, $\phi_N=\pi/3$, $\theta_L=\pi/4$, and $\phi_L=\pi/5$ in
   each case, we have assumed detection with the detector network introduced in Table 
   \ref{tab:DetCord}, and errors correspond to a fixed $\mathrm{SNR}_\mathrm{tot}=100$. 
   We find similar trends for other random choices of binary directions and orientations (not shown). }  
   \label{Table:CircLimParams}
\end{table}



\section{Summary and Conclusion} 
\label{sec:Conc} 

 We carried out a Fisher-matrix-type study to determine the accuracy with which the parameters of highly 
 eccentric BH binaries may be measured using the aLIGO-AdV-KAGRA GW detector network. Eccentricity changes 
 the GWs of binaries compared to circular binaries, in several ways. In time-domain, the gravitational 
 waveform of eccentric binaries is quasiperiodic but not sinusoidal. Relativistic precession adds a slow 
 amplitude modulation to the waveform for each polarization. Eccentricity also changes the inspiral rate 
 at which the binary separation and period shrink. We take all of these effects into account using the 
 stationary phase approximation \citep{MorenoGarridoetal1994,Mikoczietal2012}. In contrast to circular 
 binaries, the waveform of eccentric binaries includes several prominent orbital frequency harmonics, 
 general relativistic precession causes each harmonic to split into three frequencies for both GW 
 polarizations respectively, and the eccentric inspiral creates a spectrum, which is different from the
 $\tilde{h} \propto f^{-7/6}$ waveform of circular inspiral sources for each harmonic. These features in 
 the waveform make it possible to accurately determine the eccentricity and angle of periapsis, and the 
 modulated inspiral rate improves the measurement accuracy of mass parameters for eccentric inspirals. 

 The main parameters that describe eccentric inspiraling binaries are the initial pericenter distance 
 $\rho_{\mathrm{p}0}$ when the eccentricity is close to unity and the final eccentricity at the last 
 stable orbit $e_\mathrm{LSO}$. These parameters are systematically different for different formation 
 channels. Thus their measurement may have important implications on the astrophysical origin of the 
 sources
 \citep{OLearyetal2009,Cholisetal2016,Rodriguezetal2016a,Rodriguezetal2016b,Chatterjeeetal2017,Gondanetal2017,SamsingRamirezRuiz2017,SilsbeeTremaine2017}. 
 Based on a survey with $10 \, \mathrm{M}_{\odot} - 10 \, \mathrm{M}_{\odot}$ and $30 \, \mathrm{M}_{\odot} - 30 \, \mathrm{M}_{\odot}$ precessing 
 highly eccentric BH binaries at $100\, \mathrm{Mpc}$ using the planned aLIGO-AdV-KAGRA detector network,
 our results are be summarized as follows.
 \begin{enumerate}
 \item  The $\mathrm{SNR}_\mathrm{tot}$ improves by a factor of $\sim 1- 1.7$ (depending on 
 $\rho_{\mathrm{p}0}$, the component masses, and the sky position and binary orientation angles, see 
 Figure \ref{fig:SNR_Prec}) for $30 \, \mathrm{M}_{\odot} - 30 \, \mathrm{M}_{\odot}$ precessing highly eccentric BH binaries 
 compared to similar binaries in the circular limit \footnote{We adopted the leading order stationary 
 phase approximation waveform for circular sources.} with the same masses and  distance. The volume in 
 the Universe for a fixed maximum $S/N$ is \mbox{$\langle (S/N)^3\rangle \sim 5 \times$ ($2\times$)} 
 larger for eccentric inspiraling binaries with $\rho_{\mathrm{p}0}=10$ ($\rho_{\mathrm{p}0} 
 = 20$) than for similar binaries in the circular limit.
 
\item We determined how the parameters' measurement accuracies depend on the initial dimensionless
 pericenter distance ($\rho_{\mathrm{p}0}$) for precessing highly eccentric BH binaries. The smallest 
 errors are obtained for small $\rho_{\mathrm{p}0} < 10$ values for the sky position and angular 
 momentum and $\rho_{\mathrm{p}0} < 20$ for the luminosity distance and $\rho_{\mathrm{p}0}$. However, 
 the errors for fast parameters, which are sensitive to the GW phase, like the chirp mass, the initial 
 eccentricity, and the eccentricity at the last stable orbit improve most significantly for a higher 
 $\rho_{\mathrm{p}0}$ between $10$ and $80$ (Figure \ref{fig:Delta_rho}). 
 
 \item The parameter estimation errors can improve significantly for highly eccentric precessing BH binaries 
 compared to similar binaries in the circular limit by a factor of (depending on $\rho_{\mathrm{p}0}$,  
 the component masses, and the sky position and binary orientation angles, 
 see Figures \ref{fig:Dist1} and \ref{fig:Delta_rho})
\begin{itemize}
 \item  $\sim 1-200$ for the mass errors, 
 \item  $\sim 1-4.5$ for the semi-major and semi-minor axes of the sky localization ellipse, 
 \item  $\sim 1-2$ for the distance errors, 
 \item  $\sim 1-3$ for the semi-major and semi-minor axes of the error ellipse for the binary 
 orientation.
\end{itemize}

 \item \label{i:ecc} For initially highly eccentric BH binaries at $D_{\rm L}=100\,\rm Mpc$, 
 the measurement errors for parameters specific to precessing highly eccentric BH binary sources may 
 be as low as of order (depending on $\rho_{\mathrm{p}0}$, the component masses, and the sky position 
 and binary orientation angles, see Figures \ref{fig:Dist2} and \ref{fig:Delta_rho})
\begin{itemize}
 \item $10^{-5}$ for the final eccentricity errors at LSO, 
 \item $10^{-4}$ for the initial eccentricity errors,
 \item  $10^{-3}$ for the initial pericenter distance.
\end{itemize}

 \end{enumerate}

 For initially moderately eccentric to low eccentricity binaries, the parameter measurement
 errors and SNRs improve by a smaller amount for low to intermediate $\rho_{\mathrm{p}0}$ (Figure 
 \ref{fig:Err_LowEccDep}).

 Note that the eccentricity errors are remarkably low, which is not surprising given 
 that eccentricity is encoded in several measurable features of the waveform including the orbital 
 harmonics, the splitting of each frequency harmonic into triplets, the frequency evolution 
 of harmonics (the ``chirp''), the frequency evolution of the distance between the spectral 
 triplets, the low frequency cutoff of the signal at the initial pericenter frequency, and 
 the eccentricity dependence of the last stable orbit where the inspiral transitions into a
 rapid coalescence. 
 
 However, there are several factors which may significantly increase the measurement errors in
 more typical cases. First, more typical sources are expected to be at much larger distances 
 than $100\,\mathrm{Mpc}$. Assuming crudely that the eccentricity errors scale with $D_\mathrm{L}$,
 the median measurement errors for a $30 \, \mathrm{M}_{\odot} - 30 \, \mathrm{M}_{\odot}$ precessing highly eccentric BH
 binary at $\sim 410\, \mathrm{Mpc}$ (similar to GW150914) are expected to be $\Delta e_\mathrm{LSO} \sim 8.8 \times 
 10^{-4}$ ($1.3 \times 10^{-4}$) for the final eccentricity at the last stable orbit if $\rho_
 {\mathrm{p}0} = 10$ ($20$)
 for the design sensitivity of the aLIGO/AdV/KAGRA instruments. In these cases, the expected median initial eccentricity error is 
 $\Delta e_0 \sim 8.9 \times 10^{-3}$ $(1.2 \times 10^{-2})$, and the median initial pericenter 
 distance is $\Delta \rho_{\mathrm{p}0} \sim 4.4 \times 10^{-2}$ ($9.4 \times 10^{-2}$; see Table 
 \ref{tab:PrecBin}). 
 
 Another important simplifying assumption, which may have skewed the errors to lower values, was 
 to neglect higher order post-Newtonian (PN) corrections which depend on the spin of the merging
 objects. The spins of the two binary components introduce $6$ additional parameters, which may 
 become partially degenerate with all other parameters thereby increasing their errors. On the 
 other hand, spin precession breaks degeneracies between the binary orientation and other slow 
 parameters \citep{LangHughes2006,Kocsisetal2007,Chatziioannouetal2014}. However, the 
 eccentricity-induced orbital harmonics enter at the Newtonian order, the frequency-triplets 
 due to GR precession enter at the low $1$PN order. Therefore, the eccentricity-related spectral
 features are dominant already at the early stages of the inspiral when higher order PN corrections
 are negligible. For this reason, the estimated $\Delta e_0$ and $\Delta \rho_{\mathrm{p}0}$ errors
 are expected to be robust. On the other hand, most of the SNR accumulates at late times for 
 stellar-mass BH binaries where the high order PN corrections are significant. We leave an 
 estimate of parameter estimation errors for spinning binaries to future work. 
 
 Finally, an important simplification is the Fisher matrix method itself, which is valid only 
 if the waveform model is a faithful representation of the GW signal and the noise is Gaussian 
 and sufficiently small that the SNR is sufficiently large that the parameter error region is 
 an ellipsoid in parameter space and when the parameter derivate of the waveform is non-vanishing. 
 For smaller SNR, the parameter error region geometry is more complex and the uncertainties are 
 generally higher \citep[see][and reference therein]{CornishLittenberg2015}. The parameter 
 estimation errors may also be affected by theoretical uncertainties of the waveform model 
 \citep{CutlerVallisneri2007}, which may be especially important for highly eccentric binaries, 
 where the post-Newtonian expansion is slowly convergent \citep{KocsisLevin2012}. 
 
 Such low eccentricity errors may give the aLIGO-AdV-KAGRA GW detector network the capability 
 to distinguish among different astrophysical formation channels. In a companion study 
 \citep{Gondanetal2017}, we illustrate the expected distribution of eccentricities and other physical
 parameters for single-single GW capture sources in galactic nuclei. Similar studies for other 
 astrophysical formation channels are underway.

 Future multi-waveband searches of eccentric inspiral sources with \textit{LISA} and aLIGO-AdV-KAGRA
 \citep{KocsisLevin2012,Sesana2016} have good prospects for even more accurate measurements of the
 physical parameters well beyond the level reported here. In that case the GW frequency range is 
 much wider. Since eccentricity decreases due to GW emission, eccentricity may be expected to be 
 much higher at lower frequencies in the \textit{LISA} band. This leads to a much larger total GW 
 phase shift caused by eccentricity. The better measurement of relativistic precession may more 
 efficiently break degeneracies between mass and other parameters. The modulation caused by the 
 orbit of the instrument around the Sun and Earth's spin can help break degeneracies among source
 direction, orientation, and other parameters. Accounting for eccentricity for third generation 
 Earth based (e.g. Einstein telescope), deci-Hertz to mHz space-based instruments  will be essential 
 \citep{ChenAmaroSeoane17}.
 
 
 
\section*{Acknowledgment}
 We thank the anonymous referee for constructive comments, which helped improve the quality of the 
 paper. This project has received funding from the European Research Council (ERC) under the European
 Union's Horizon 2020 research and innovation programme under grant agreement No 638435 (GalNUC) and 
 by the Hungarian National Research, Development, and Innovation Office grant NKFIH KH-125675. This 
 work was performed in part at the Aspen Center for Physics, which is supported by National Science 
 Foundation grant PHY-1607761. The calculations were carried out on the NIIF HPC cluster at the
 University of Debrecen, Hungary.
 
 
 
