text:
\section{Introduction}
\label{Sect:intro}

The $\phi$ meson is a particularly interesting hadron for a variety of reasons.
It is a non-strange particle composed of a $\text{s}\bar{\text{s}}$ pair of quarks,
and is characterized by the narrow mass distribution
centered at $m_{\phi} = 1.0195$~GeV/$c^2$, and decay width
$\Gamma = 4.27$~MeV~\cite{PDG}.
Its suppressed decay modes have accompanied the discovery of the
Okubo-Zweig-Iizuka (OZI) rule~\cite{OZI}.
An enhancement of its production in the ultra-relativistic heavy-ion (AA)
collisions has been interpreted as the result of fragmentation of gluons
into q$\bar{\text{q}}$ pairs in the Quark-Gluon Plasma~\cite{Koch90}.
Calculations in the frame of the transport BUU code suggest that $\phi$ mesons
from the heavy-ion collisions at energies below threshold in a free
nucleon-nucleon (NN) collision may be produced on average even earlier
than kaons, and due to relatively weak absorption and rescattering
cross sections, a majority of $\phi$'s may survive the collision.
Therefore, $\phi$ mesons may be a good probe of earlier stages of
the heavy-ion collision~\cite{Scha10}.

In the proton-proton (pp) collisions, the ratio of $\phi$ to $\omega$
cross sections obtained near $\phi$ meson threshold was found
to violate the OZI rule by about an order of magnitude~\cite{Bale01,Hart06}.
Studies of near-threshold $\phi$ meson production in pA collisions
revealed the target mass dependence of A$^{0.56(3)}$~\cite{Taba06},
and pointed to the momentum-dependent widening of the $\phi$ width
in the nuclear medium, from where the $\phi$N absorption cross
section of about 15--20~mb was estimated~\cite{Hart12}.

In the domain of heavy-ion collisions data on the subthreshold $\phi$ meson
production is quite scarce. In the first measurement only $23 \pm 7$
events were found in the K$^+$K$^-$ channel
and the reconstructed yield per event varied by the \mbox{factor 4}
depending on the value of temperature parameter of the $\phi$
meson source assumed for the efficiency evaluation~\cite{Mang03}.
Recently, three samples of about \mbox{110--170} $\phi$ mesons were measured
at quite similar values of the average number of participant nucleons
($\langle A_{\text{part}} \rangle_{\text{b}}$ = \mbox{35--50}),
and beam kinetic energies in the range of 1.75--1.9A GeV
~\cite{Agak09,Lore10,Pias15,Gasi16}. Production yields per triggered event
were found to be $2.5 \text{--} 4.5 \times 10^{-4}$.
Taking into account that the branching ratio equals
$\text{BR} (\phi \rightarrow \text{K}^+\text{K}^-$) = 48.9 \%~\cite{PDG},
and comparing to the corresponding production yields of K$^-$,
it was inferred that $\phi$ mesons were the source of about
\mbox{15--20~\%} of negatively charged kaons.
A determination of the size of this contribution,
and tracing its dependency on the beam energy and collision centrality
is important in the context of the investigation of in-medium modifications
of K$^-$ mesons from their emission patterns~\cite{Fuch06,Lutz04,Hart11,Wisn00,Fors07}.
As the mean decay path of $\phi$ meson is 46~fm, it decays
mostly outside the collision zone, and therefore the majority of
daughter kaons are not produced in the medium.
Phase space distribution of negative kaons is thus composed
of two contributions, and the extraction procedure of
the kaon potential inside the medium should take into account
the side feeding from $\phi$ meson decays.

In recent years there has been a considerable development
on the side of transport models.
Earlier calculations in the frame of the BUU code
performed for the collisions of central Ni+Ni at 1.93A GeV,
and Ru+Ru at 1.69A GeV
favored the dominance of meson-baryon production channels
($\text{MB} \rightarrow \phi \text{B}$, M = [$\rho$,~$\pi$],
B = [N,~$\Delta$])~\cite{Barz02}.
However, further BUU calculations employing the new effective
parameterization of the
$\text{pN} \rightarrow \text{pN} \phi$ channels (N = [p,n]),
performed for the central Ni (1.93A GeV) + Ni and semi-central
Ar (1.756A GeV) + KCl collisions
suggested the BB production channels to be of equal strength as MB.
An interesting study of $\phi$ meson production in Ca (1.76A GeV) + Ca
collisions in the frame of the UrQMD model has been recently published~\cite{Stei16}.
Within this approach the decays of massive resonances
[$N^*(1990) ~...~ N^*(2250)$] have been pointed out as a
dominant source of $\phi$ mesons (cf. Fig.~6 therein).
In both types of calculations the \mbox{$\phi\slash\text{K}^-$}
ratio was predicted to be similar to the experimental findings
in the Ar+KCl and Ni+Ni systems.
The BUU model gave detailed predictions on the composition of production channels
as a function of impact parameter, distribution of production times
and kinematic variables, while the UrQMD calculations predicted a specific
profile of the excitation function of \mbox{$\phi\slash\text{K}^-$} ratio
with a maximum around the threshold energy,
evolving around the AGS beam energy region into the constant, lower value.
Thermal models also deliver predictions of the $\phi$ meson
production yields in an event.
In Ref.~\cite{Cley99} the \mbox{$\phi\slash\pi^+$} ratio is predicted
to be centrality independent,
and correlated only with the temperature of the system.
Calculations of the \mbox{$\phi\slash\text{K}^-$} ratio shown in Fig.~13 of
Ref.~\cite{Agak09} predict the strong enhancement around the threshold,
sensitive to the suppression of the volume in which the open strangeness
is produced ($R_{\text{C}}$ parameter).

However, there seems to be a strong disparity between the variety
of predictions of subthreshold $\phi$ meson production
and the scarcity of available experimental data.

In our work we address this gap by analysing $\phi$ mesons produced
in the central Ni+Ni collisions at the beam kinetic energy of 1.93A GeV,
and decaying in the dominant K$^+$K$^-$ channel.
While, as mentioned above, these mesons were already investigated
within the similar centrality class,
a relatively low number of $4.7 \times 10^6$ acquired events
in the experiment reported in~\cite{Mang03}
resulted in only $23 \pm 7$ events attributed to $\phi$ mesons.
In the experiment described in this paper, about 17 times more events
were acquired for the central trigger (defined in Sect.~\ref{Sect:exp}),
giving a good chance to measure $\phi$ mesons at considerably
better significance level.

As data on $\phi$ meson production at 1.9A GeV was also published
for the less central Ni+Ni collisions, as well as central Al+Al collisions,
our investigation allows for the first time to gain insight into
the centrality dependence of the subthreshold $\phi$ meson production
in heavy-ion collisions.

\section{Experiment}
\label{Sect:exp}

The S261 experiment was performed at the SIS--18 synchrotron in the GSI, Darmstadt.
$^{58}$Ni ions were accelerated to the kinetic energy of 1.93A GeV,
and were incident on the $^{58}$Ni target of the FOPI modular spectrometer
with the average intensity of $4\text{--}5 \times 10^5$ ions per spill.
The target had thickness of 360 mg/cm$^2$
corresponding to 1.0\% interaction probability.

The innermost detector of the FOPI setup is the
Central Drift Chamber (CDC) covering the wide range
of polar angles ($27^\circ < \theta_{\text{lab}} < 113^\circ$,
\footnote{The angles are given with respect to the target position
in the S261 experiment.}),
and subdivided azimuthally into 16 identical sectors.
The CDC was surrounded by the Time-of-Flight (ToF) detector named
Plastic Scintillation Barrel (PSB), which
covered the polar angles of $26.5^\circ < \theta_{\text{lab}} < 56^\circ$.
Both detectors were mounted inside the magnet solenoid,
which generated the field of 0.6~T. The Plastic scintillation
Wall (PlaWa) was placed downstream of the CDC.
More information on the geometry and performance of the FOPI apparatus
can be found in~\cite{FOPI}.

The on-line "central" trigger was based on the large number of hits
in the PlaWa.
After an additional off-line removal of events with vertex position
outside the target, about $79\times 10^6$ events were selected,
amounting to $(22 \pm 1)$\% of the total reaction cross-section.
Assuming the sharp cut-off approximation between
the total reaction cross section and maximum impact parameter,
and the geometrical model of interpenetrating spheres,
the number of participant nucleons averaged over the impact parameter
was estimated to be $\langle A_{\text{part}} \rangle_{\text{b}} = 74 \pm 2$.

\section{Data analysis}
\label{Sect:ana}

Charged particles passing through the CDC detector
activate the nearest sense wires,
represented within the data acquisition scheme as "hits".
An off-line tracking routine collects them into "tracks",
and reconstructs the emission angles and specific energy loss of particles.
For these particles which additionally generate hits in the PSB detector,
the ToF information is also obtained.
A reconstruction of the collision vertex allows for the elimination of
events not originating from the target.

To preselect the tracks with good reconstruction quality,
these composed of few hits were rejected.
The tracks with too large distance of closest
approach to the vertex were also eliminated.
The candidates for K$^+$ and K$^-$ mesons were searched
from the tracks matched with hits in the ToF barrel.
For a good quality of matching, the extrapolation of the
CDC track was required to pass close to the hit candidate
in the PSB detector.
To minimize the edge effects, the range of accepted polar angles
was trimmed down to $30^\circ < \theta_{\text{lab}} < 53^\circ$,
and K$^+$ (K$^-$) candidates were required to have the transverse
momentum $p_{\text{T}}$ of at least 0.18 (0.14) GeV/$c$.
A motivation for these cuts will be also discussed in Sect.~\ref{Sect:eff}.

\begin{figure}[tbh]
 \includegraphics[width=8.6cm]{s261_mcdc-mpsb.eps}
 \caption{\label{fig:kpkmmass}(Color online)
  Identification of charged kaons on the
  two-dimensional map of FOPI mass parameters: $m_{\text{CDC}}$
  (obtained from the CDC) versus $m_{\text{PSB}}$
  (from the combination of CDC-PSB data). Rectangles in dashed lines
  delineate the cuts imposed for the selection of K$^+$ and K$^-$ candidates.}
\end{figure}

Within the basic identification technique of charged particles in
FOPI the information from a drift chamber (here: the CDC) is used.
A sign of electric charge is inferred from the direction
of the track curvature on the transverse plane.
Subsequently, the momentum and specific energy loss inside the chamber
is substituted to the Bethe-Bloch formula, and the mass parameter
is extracted (dubbed here $m_{\text{CDC}}$).
However, the resolution of the specific energy loss is too weak for
the identification of kaons in a reasonable momentum range.
Therefore, an additional identification technique is applied:
the mass is reconstructed by substituting the momentum and velocity
to the relativistic formula $p = m \gamma v$. The mass parameter
obtained in this method is dubbed $m_{\text{PSB}}$.
For the present analysis both methods were applied together.
In order to reject most of the pion and baryon events, the gates
were imposed on the abovementioned mass parameters,
shown as dashed rectangles in Fig.~\ref{fig:kpkmmass}.

The $\phi$ mesons in the K$^+$K$^-$ decay channel were investigated
in the invariant mass ($M_{\text{inv}}$) distribution of kaon pairs.
As the resolutions of $m_{\text{CDC}}$ and $m_{\text{PSB}}$ distributions
in FOPI are well known to deteriorate with momentum,
to limit the influence of particles with neighbouring masses
(pions and protons), the maximum momentum was set to 0.65~GeV/$c$.
In the FOPI analyses aiming at the selection of K$^+$ and/or K$^-$
using PSB, a reasonable level of signal to background
(S/B) ratio was usually kept if particles were limited to lower
maximum momenta: about 0.55~GeV/$c$ for K$^+$,
and 0.38~GeV/$c$ for K$^-$~\cite{Wisn00,Ziny14,Pias15,Gasi16}.
However, for the $\phi$ meson reconstruction analysis an unambiguous
selection of kaons is not so critical, as some admixture of misidentified
particles is not expected to create a correlation around $m_{\phi}$.
To minimize possible distortions arising for the kaon pair candidates with
small relative angles, pairs were required not to intersect in the active
region of the CDC.
The background of uncorrelated pairs was reconstructed using
the mixed events method, where K$^+$ and K$^-$ tracks from different
events were matched. The events from where two kaons were paired
had to be attributed to the same centrality class,
determined by the multiplicity of tracks in the CDC.
To prevent from smearing out of possible flow pattern of kaons
by matching events with randomly oriented reaction planes,
the parent events were rotated azimuthally to align these planes.
Finally, the background was normalized to the true pair distribution
in the region $1.05 < M_{\text{inv}} < 1.18$~GeV/$c^2$.

\begin{figure}[tbh]
 \includegraphics[width=8.6cm]{s261_phi_minv.eps}
 \caption{\label{fig:phiminv}(Color online)
   (a) Invariant mass plot of true (solid line),
   and mixed (scattered points) K$^+$K$^-$ pairs. The normalisation window
   for the spectrum of mixed pairs is indicated with the arrow.
   (b) $\phi$ meson signal obtained after the background subtraction.
   The solid curve shows the Gaussian distribution fitted within
   $1.00 < M_{\text{inv}} < 1.04$~GeV/$c^2$.
   (c) A comparison of profiles of experimental (points) and simulated
   (solid line) $\phi$ meson signals.}
\end{figure}

The distribution of true and mixed pairs is shown in
Fig.~\ref{fig:phiminv}(a). The $\phi$ meson signal obtained
by the subtraction of background from the spectrum of true pairs is
presented in panel (b) and exhibits some excess of counts in the region
$M_{\text{inv}} < 1.01$~GeV/$c^2$. To check if this excess
is generated by true $\phi$ meson events or if it is a side effect,
the experimental signal was compared to the profile obtained
in the simulation of $\phi$ meson production and detection,
shown by the solid line in Fig.~\ref{fig:phiminv}(c).
Here we only report the final result of this simulation,
as it will be described in details further in Sect.~\ref{Sect:eff}.
As both distributions exhibit the abovementioned excess,
the $\phi$ meson signal was interpreted as the total
number of signal counts in the range $0.98 < M_{\text{inv}} < 1.05$~GeV/$c^2$.
Within the set of cuts described above,
the number of identified $\phi$ mesons was found to be $110 \pm 19$.
In addition, a Gaussian distribution was fitted to the signal
in the range of $1.00 < M_{\text{inv}} < 1.04$~GeV/$c^2$,
yielding the $\sigma$ parameter of 6~MeV/$c^2$
(see Fig.~\ref{fig:phiminv}(b)).
The peak maximum was found to be $1.0211(13)$~GeV/$c^2$,
about 1 standard deviation away from the nominal value in vacuum.
Within the range of $\pm ~ 2\sigma$ the signal to background ratio
was found to be 0.4, and the significance 5.0.
The number of reconstructed $\phi$ mesons and the other
characteristics of $M_{\text{inv}}$ distribution depend on the choice
of cuts applied to the track quality parameters, mass parameters,
maximum momenta, binning, and windows for integration and normalization.
This issue is further discussed in Appendix~\ref{Sect:detaSystErr}.

\section{Efficiency determination}
\label{Sect:eff}

The calculation of efficiency was performed using the GEANT3-based
environment~\cite{GEANT}, where the modeled detectors were
positioned as in the S261 experiment. \mbox{$\phi$ mesons} were sampled
from the relativistic Breit-Wigner mass profile with the nominal mass
and decay constant~\cite{PDG}, and the phase space was populated
according to the Boltzmann function multiplied by the
simple anisotropy term:
\begin{equation}
 \label{EQ:phigeantmodel}
 \frac{d^2N}{dEd\theta} \sim ~pE~ \exp (-E\slash T_{\text{s}})  \cdot \left(1 + a_2 \cos^2 \theta\right) ,
\end{equation}
\noindent where $T_{\text{s}}$ is the temperature of the source,
and $a_2$ is the anisotropy parameter. Little is known about
the characteristics of the $\phi$ meson distribution in phase space.
The inverse slopes for the collisions at 1.9A GeV, and lower
$\langle A_{\text{part}} \rangle_{\text{b}}$, were found to be
between 70 and 130 MeV~\cite{Pias15,Gasi16}.
As for baryons and K$^\pm$ mesons more central collisions translate
into somewhat higher inverse slopes~\cite{Hong97,Fors07},
one may cautiously apply this trend to $\phi$ mesons.
Thus, for the determination of the efficiency in the present analysis
$T_{\text{s}}$ was varied between 80 and 140 MeV.
Regarding the polar anisotropy of $\phi$ meson emission,
no data is known so far.
Therefore, as roughly expected from the systematics for other strange
particles produced in the studied beam energy range and collision
centrality~\cite{Fors07}, $a_2$ was varied between 0 and 1.
The generated $\phi$ mesons were boosted to the laboratory frame,
and all of them were allowed to decay into K$^+$K$^-$ pairs
(the BR factor was accounted for in the calculation of the efficiency).
Particles were subsequently added to the events
of Ni+Ni collisions generated by the IQMD transport code~\cite{IQMD},
which aims at providing the realistic set of particles emitted
from the heavy-ion collisions. However, kaons emitted directly from the
collision zone were not produced by this code.

After the propagation and digitization stages performed by the GEANT routines,
events were processed using the same tracking, matching and
correlation routines, as for the experimental data.
Distributions of the quality parameters for the experimental and simulated
tracks were found to mostly overlap. To account for minor discrepancies,
the cut values for the simulated tracks were tuned such that
a given cut filtered out the same fraction of tracks in the
generated sample, and in the true (measured) one.
However, the distributions of mass parameters for the
tracks in the simulation were found to be narrower
than these from the experiment.
To tackle this problem, in the first step
for the experimental tracks of K$^+$ and K$^-$ candidates
the two-dimensional distributions of $m_{\text{CDC}}$ versus $m_{\text{PSB}}$
were inspected for consecutive slices of momentum,
as shown in Fig.~\ref{fig:kpkmmass}.
In the next step, two-dimensional Gaussian functions were fitted
to these distributions within regions not influenced
by other particles (pions and/or protons).
A comparison of the extracted Gaussian distribution
to the rectangular cuts applied on the $m_{\text{CDC}}-m_{\text{PSB}}$
plane (see dashed rectangles in the abovementioned Figure)
allowed to find the fraction $f(p)$ of kaons
rejected by these cuts. Depending on momentum and particle sign,
$f(p)$ ranged from 0.01 to 0.15.
On the side of the simulation, the narrower mass spectra allowed to
select the sample of kaons nearly unambiguously. However, to reproduce
the abovementioned loss of kaons occuring for the experimental data,
a fraction $f(p)$ of simulated kaons was removed from the sample.

The obtained profile of the invariant mass spectrum of K$^+$K$^-$
pairs is shown as solid line in Fig.~\ref{fig:phiminv}(c).
The \mbox{$\phi$ meson} reconstruction efficiency in the K$^+$K$^-$ decay channel,
calculated for the set of cuts described above, and for the source parameters
$T_{\text{s}} = 120$~MeV, and $a_2 = 0$ was found to be $3.5 \times 10^{-3}$.
It was also found that the efficiency obtained with $a_2 = 1.0$
was only 2\% lower, while variations due to different assumed values of
$T_{\text{s}}$ remained within $\pm 10 \%$ range.

In the next step we have checked whether the simulation reproduced
the efficiency of matching of tracks in the CDC with hits in the PSB.
As the statistics on kaons was too low for a reliable
study of this issue, the effect for negatively charged particles
was tested on $\pi^-$ tracks, and that for positively charged ones
-- on the sum of $\pi^+$ and protons.
In the first step the ratio of CDC tracks with associated
hit in the PSB to all the reconstructed CDC tracks was inspected
as a function of laboratory polar angle and transverse momentum.
This ratio was obtained independently for the experimental
and simulated data. Subsequently, the ratio obtained for
the experimental data was divided by that for the simulated data,
to yield the correction factor for the matching efficiency,

\begin{equation}
\varepsilon^{CDC-PSB} \left( \theta_{\text{lab}}, p_{\text{t}} \right)
 = \frac{N^{\text{PSB}}_{\text{exp}}}{N^{\text{CDC}}_{\text{exp}}}
   \slash
   \frac{N^{\text{PSB}}_{\text{sim}}}{N^{\text{CDC}}_{\text{sim}}} \quad .
\end{equation}

\noindent The meaning of this factor is:
if the matching efficiency for a particle is worse (better) in the experiment
than in the simulation, the factor is lower (higher) than 1.
The $\theta_{\text{lab}} \text{--} p_{\text{T}}$ distribution of
$\varepsilon^{CDC-PSB}$,
presented in Fig.~\ref{fig:bareff}, shows that discrepancies between the
simulated and experimental matching efficiencies remain
mostly within a $\pm 8\%$ range.
However, they arise at the detector edges and at low transverse momenta.
This finding was a motivation to limit the range of accepted
polar angles to $30^\circ < \theta_{\text{lab}} < 53^\circ$,
and set the minimum $p_{\text{T}}$ at 0.18 (0.14) GeV/$c$ for K$^+$ (K$^-$),
as mentioned above in Sect.~\ref{Sect:ana}.
Still, to correct for the remaining discrepancies,
the kaon candidates selected to build the experimental
$M_{\text{inv}}$ spectrum were weighted with the factor
\mbox{$1 \slash \varepsilon^{CDC-PSB} (\theta_{\text{lab}}, p_{\text{t}})$}.

More details on the estimation of systematic errors, conventions
applied and correction of the $\langle A_{\text{part}} \rangle_{\text{b}}$
value for the experiment presented in~\cite{Pias15} is reported in
Appendix \ref{Sect:deta}.

\begin{figure}[bt]
 \includegraphics[width=8.6cm]{s261_bareff_pim-ppip.eps}
 \caption{\label{fig:bareff}Maps of correction factor for possible
unequal efficiencies of matching of CDC tracks with PSB hits between
the experimental and simulated data,
obtained for (a) $\pi^-$ tracks, and (b) sum of $\pi^+$ and proton tracks.
Regions in between dashed lines correspond to cuts applied in
the K$^+$K$^-$ analysis.}
\end{figure}

\section{Results}
\label{Sect:resu}

\begin{table}
 \caption{Available data on the centrality dependence of the $\phi$ meson
 production yield per triggered event at 1.9A GeV.
 Note, that some of published results were slightly adjusted
 according to the conventions defined in Appendix \ref{Sect:deta}.}
 \label{Tab:PhiYields}
 \begin{tabular}{cccc}
 \hline
 Reaction & $\langle A_{\text{part}} \rangle_{\text{b}}$ & Yield $\times 10^{-4}$ & Ref. \\
 \hline
  Al+Al   &   42.5 & $3.1 \pm 0.5 \text{(stat)} \pm 0.3 (\text{syst})$ & \cite{Gasi16} \\
  Ni+Ni   &   46.5 & $4.4 \pm 0.7 \text{(stat)} \pm 1.3 (\text{syst})$ & \cite{Pias15} \\
  Ni+Ni   &   74   & $8.6 \pm 1.6 \text{(stat)} \pm 1.5 (\text{syst})$ & This work \\
 \hline
  Ni+Ni (BUU) & 86 & 25. & \cite{Scha10} \\
 \hline
 \end{tabular}
\end{table}

\begin{figure}[b]
 \includegraphics[width=8.4cm]{phi_Apart_fit.eps}
 \caption{\label{fig:phiApart}(Color online)
  $\phi$ meson yield per event at the beam energy of 1.9A GeV
  as a function of number of participant nucleons averaged
  over the impact parameters.
  Full circles: data measured at 1.9A GeV;
  open circle: data obtained at 1.756A GeV~\cite{Agak09}
  and scaled according to the systematics in~\cite{Meta93};
  Dotted curve: power function fit
  to the points marked by full circles. See text for details.
 }
\end{figure}

The production yield of $\phi$ mesons per triggered
collision of Ni+Ni ions was found to be:
\begin{equation}
 \rm{P(\phi) = [8.6 ~\pm~ 1.6 ~(\text{stat}) \pm 1.5 ~(\text{syst})]
     \times 10^{-4}}\quad .
 \label{eq:yield_phi}
\end{equation}

\noindent where the systematic error is given within the 68.3\%
confidence level (see Appendix~\ref{Sect:detaSystErr} for details).
For the consistency check we have verified that if the same
centrality cut is imposed on the events of Ni+Ni collisions
reported in~\cite{Pias15}, the $\phi$ meson yield extracted
from that data agrees within 1 standard deviation with the
abovementioned result.

A compilation of the yields per event obtained by FOPI at 1.9A GeV
is shown in Tab.~\ref{Tab:PhiYields}
and in Fig.~\ref{fig:phiApart} (full circles).
The multiplicity of $\phi$ mesons from Ar+KCl collisions at 1.756A GeV,
obtained by the HADES collaboration,
cannot be compared directly to this sample,
as due to the scarcity of data on subthreshold $\phi$ meson production
the excitation function is not known.
However, a common systematics of the excitation function of
$\pi^{0,-}$, K$^+$, $\eta$, and $\rho^0$ meson yields per participant
nucleon has been found, as shown in Fig.~6 of Ref.~\cite{Meta93}.
Interestingly, if the $\phi$ meson yield from Ar+KCl collisions
is scaled using this systematics, the ''adjusted'' value
shown as open circle in Fig.~\ref{fig:phiApart} falls well
in line with the FOPI data points obtained at 1.9A GeV.

We have fitted the FOPI data with the power function,
\mbox{$P (\phi) \sim A_{\text{part}}^\alpha$},
and obtained $\alpha = 1.8 \pm 0.6$.
For the purpose of fitting, the statistical and systematic errors
were added quadratically.
Within the thermal model approach without the strangeness undersaturation
effects the volume dependence of the $\phi$ meson multiplicity in an event
is predicted to be linear~\cite{Cley99}.
The measured $\alpha$ parameter is in agreement with this prediction
within 1.4 standard deviations.
According to the calculations in the frame of the BUU transport code
the $\phi$ meson yield in an event should be equal to $2.5 \times 10^{-3}$
for the class of 12\% most central Ni+Ni collisions
(corresponding to $\langle A_{\text{part}} \rangle_{\text{b}} = 86$
within the sharp cut-off approximation)~\cite{Scha10}.
The experimental data for this centrality class (see~\cite{Mang03})
is too scarce for a precise comparison with this prediction.
However, the BUU result appears to overpredict the abovementioned
power function, if the latter is slightly extrapolated beyond the region
of experimental data points.

The obtained centrality dependence of the $\phi$ meson production yield
can be compared to the $A_{\text{part}}$ scaling of the other mesons
produced at similar beam energies.
The pattern for $\pi^{+,-}$ mesons exhibits somewhat
sublinear dependence~\cite{Pelt97,Reis07}.
For K$^0_{\text{s}}$, $\alpha = 1.20 \pm 0.25$ was reported~\cite{Mers07}.
Turning to the charged kaons, the centrality dependence of their
multiplicity was not measured at 1.9A GeV,
but at 1.8A GeV the exponent was found to be
$\alpha = 1.65 \pm 0.15$ for K$^+$, and
$\alpha = 1.8 \pm 0.3$ for K$^-$~\cite{Bart97}.

These observations can be interpreted by a mixture of different mechanisms.
At beam energies around 1.9A GeV charged pions are produced
abundantly far above the threshold in the free NN collision,
with considerable feeding from $\Delta$ resonance decays~\cite{Hong97}.
Their sublinear centrality dependence was interpreted
as the result of interplay between the surface ($\sim \text{A}^{2/3}$)
and volume ($\sim$A) contributions~\cite{Reis07}.
On the other hand the analysis of $\pi^+$ mesons emitted from
the collisions around \mbox{0.8--1.0A~GeV}
has shown that whereas for pions with energies allowed in the
free NN channel the $\alpha$ exponent is consistent \mbox{with 1},
for high energy pions it exceeds the unity. In the latter domain
the exponent rises with the pion energy,
as more and more additional energy is necessary to pass the
threshold in the NN channel (see Fig.~11 in Ref.~\cite{Muen97}).

The rise of $\alpha$ coefficient with probing more subthreshold energies
was also observed for the K$^+$ production~\cite{Bart97,Fors07}.
At the beam energies of 1--2A GeV these particles are produced
around the threshold in the $\text{NN} \rightarrow \text{NKY}$ channel,
which is energetically the lowest-lying channel
for the K$^+$ production in the free NN collision.
Therefore the multi-step channels involving the intermediate
$\Delta$ baryons or pions, which usually facilitate the subthreshold
production of particles, are the important competitors,
and -- according to the BUU and IQMD transport codes --
these channels are responsible for
over a half of the kaon production cross section~\cite{Scha10,Hart11}.
Interestingly, within the thermal model approach the strong rising trend
of K$^+$ yield in an event with number of participants
is explained by the associated production of a hyperon in the
$\text{NN} \rightarrow \text{NKY}$ channel, required
to conserve the strangeness~\cite{Cley99}.
Regarding the K$^-$ production, at beam energies between 1 and 2A GeV
they are produced quite deeply below the threshold in the free NN
collision, and the centrality dependence of K$^-$ multiplicity in an event
was found to be stronger than linear, repeating the case of K$^+$ and pions.
In addition, for a fixed energy the centrality dependence of K$^-$ yield
appeared to follow the pattern for K$^+$~\cite{Bart97,Fors07}.
The latter observation can be explained by noticing that the
strangeness-exchange channels ($\pi\text{Y} \rightarrow \text{BK}^-$ and
$\text{BY} \rightarrow \text{BBK}^-$),
which are expected to dominate the K$^-$ production
(cf.~\cite{Fuch06,Scha10,Hart11}), require the presence of a hyperon, while
the energetically lowest-lying channel leading to the hyperon production,
$\text{NN} \rightarrow \text{NYK}^{+,0}$, generates also the kaon.

Turning back to the $\phi$ meson production, the extracted value
of the exponent is characterized by the large uncertainty.
Therefore, at the present level of statistics, the production mechanism
cannot be judged from the comparison of $\alpha$ parameters
between $\phi$ and other mesons.
Interestingly, the state-of-the-art calculations in the
frame of the BUU and UrQMD models advocate the entirely different
mechanisms of the $\phi$ meson production~\cite{Scha10,Stei16}.

\begin{figure}[b]
 \includegraphics[width=8.6cm]{phikm_phipip.eps}
 \caption{\label{fig:phikmpip}(Color online)
  (a) Centrality dependence of \mbox{$\phi\slash\text{K}^-$} ratio
      at the beam energy of 1.9A GeV.
      The horizontal line and the arrow mark the average value.
  (b) Centrality dependence of \mbox{$\phi\slash\pi^+$} ratio.
      The dotted line shows the fit of the constant function, and
      the dashed line -- the power function fit.
      The value of $\alpha$ exponent is the result of the latter fit.
      See text for details.
}
\end{figure}

Considering the ratio of \mbox{$\phi\slash\text{K}^-$} yields,
in case of the central Ni+Ni collisions,
we notice that the sample of negative kaons reported in~\cite{Menz00}
was measured within nearly the same centrality class as
$\phi$ mesons from the S261 experiment.
To obtain the production yield per event of these kaons,
we have integrated the rapidity distribution shown in Fig.~2 of
Ref.~\cite{Menz00}. It allowed us to find the ratio of yields:

\begin{equation}
 \rm{\frac{P(\phi)}{P(\text{K}^-)} = 0.41 \pm 0.08 (stat) \pm 0.08 (syst)
     \quad ,}
\end{equation}

\noindent where the systematic error is given within the 68.3\%
confidence level and includes the relevant uncertainties associated
to both the experiments. With this information
we have plotted the centrality dependence of the
\mbox{$\phi\slash\text{K}^-$} ratio, as shown in Fig.~\ref{fig:phikmpip}(a).
No clear trend is observed. The average value of this ratio is
$0.36 \pm 0.05$. Taking into account the branching ratio of
the $\phi \rightarrow \text{K}^+\text{K}^-$ decay we conclude
that at the beam energy of 1.9A GeV,
$(18 \pm 3)$ \% of negative kaons originate from the decays
of $\phi$ mesons, in agreement with the previous analyses~\cite{Pias15,Gasi16}.
Interestingly, these results are the same as for the Ar+KCl case,
measured at the more subthreshold beam energy of 1.756A GeV~\cite{Agak09}.
Within the version of statistical model in which the canonical ensemble
is applied to particles with open strangeness, the studied ratio is
sensitive to the correlation length $R_{\text{C}}$,
which corresponds to the reduced volume where strangeness can be produced.
As shown in Fig.~13 of Ref.~\cite{Agak09}, this sensitivity becomes
particularly strong around the $\phi$ and K$^-$ threshold energies.
The value of \mbox{$\phi\slash\text{K}^-$}
ratio obtained in the present analysis, as well as from the abovementioned
Ar+KCl collisions fall within the range of $R_{\text{C}} \in (2.2, 3.2)$~fm.
Most recent calculations within the UrQMD model appear to reproduce the found
\mbox{$\phi\slash\text{K}^-$} ratio as well (cf. fig.~4a in~\cite{Stei16}).

Fig~\ref{fig:phikmpip}(b) shows the centrality dependence of the
\mbox{$\phi\slash\pi^+$} ratio (see Appendix~\ref{Sect:detaCenKmPi}
for more details).
We have fitted the data with the constant function
(shown as dotted line) and obtained the value
$(8.3 \pm 1.1) \times 10^{-5}$ at \mbox{$\chi^2\slash\text{NDF} = 1.3$}.
An alternative fit of the power function,
\mbox{$P (\phi) \slash P(\pi^+) \sim A_{\text{part}}^\alpha$},
yields $\alpha = 1.0 \pm 0.6 ^{+0.1}_{-0.0}$ at
\mbox{$\chi^2\slash\text{NDF}$}~=~0.3.
A minor systematic error accounts for some uncertainty of the
$\pi^+$ yield at $\langle A_{\text{part}} \rangle_{\text{b}} = 42.5$,
which is discussed in Appendix~\ref{Sect:detaCenKmPi}.

The statistical model predicts the \mbox{$\phi\slash\pi^+$} ratio to be
independent of the number of participant nucleons
(see Fig.~2 in~\cite{Cley99}).
On the other hand, according to the most recent calculations within
the UrQMD model, the predicted trend of the abovementioned ratio,
shown in fig.~4b in Ref.~\cite{Stei16}, is strongly rising with $A_{\text{part}}$.
Judging by the obtained values of $\chi^2\slash\text{NDF}$, both
scenarios are consistent with our data.
Within the latter model the \mbox{$\phi\slash\pi^+$}
ratio is predicted to be sensitive to the temperature of the collision zone,
which by comparison to the abovementioned figure in~\cite{Cley99},
would be estimated to be around 70~MeV.

\section{Summary}
We have investigated the $\phi$ meson production in the K$^+$K$^-$
decay channel in the 22\% most central collisions of Ni+Ni
at the beam kinetic energy of 1.93A GeV
and found the sample of $110 \pm 19$ mesons.
After the efficiency corrections the yield per triggered event was found to be
$[8.6 ~\pm~ 1.6 ~(\text{stat}) \pm 1.5 ~(\text{syst})] \times 10^{-4}$.
It allowed for the first time to gain insight into the centrality
dependence of the multiplicity of $\phi$ mesons produced
in the heavy-ion collisions below the NN threshold.
Its dependence can be described by a power law with the exponent
$\alpha$ of \mbox{$1.8 \pm 0.6$}.
The current BUU calculations seem to overpredict this trend.
The \mbox{$\phi\slash\text{K}^-$} ratio as a function of $A_{\text{part}}$
does not exhibit any trend within the available accuracy,
and the average value was found to be $0.36 \pm 0.05$.
We have also found that the \mbox{$\phi\slash\pi^+$} ratio
has an average value of $(8.3 \pm 1.1) \times 10^{-5}$.
A power function parameterization of the $A_{\text{part}}$ dependence
of this ratio yielded the exponent of $\alpha = 1.0 \pm 0.6 ^{+0.1}_{-0.0}$.
According to the values of $\chi^2\slash\text{NDF}$ for these fits,
both of them are consistent with the data.

\begin{table}
 \caption{Variations of the cut parameters used in the investigation of
          $\phi$ mesons in the central Ni (1.9A GeV) + Ni collisions}
 \label{Tab:Syst}
 \begin{ruledtabular}
  \begin{tabular}{p{1.7cm}|p{4.6cm}|p{3cm}}

Parameter
          & Meaning
          & Values \\
\hline
$T$ [MeV] & Temperature of $\phi$ meson source
          & 80,~100, ...,~140 \\
$a_2$     & Anisotropy parameter of $\phi$ meson source
          & {0.0 , 1.0} \\
$M_{\text{inv,norm}}^{\text{min}}$ [GeV/$c^2$]
          & Lower endpoint of range of normalisation of background to true pairs
          & 1.05, 1.10 \\
$M_{\text{inv,int}}^{[\text{min--max}]}$ [GeV/$c^2$]
          & Integration window of $\phi$ meson signal in the $M_{\text{inv}}$ spectrum
          & \mbox{[0.98~--~1.05]}, \mbox{[1.00~--~1.04]}, \mbox{[1.01~--~1.04]} \\
$N_{\text{bins}}$
          & No. of bins for $M_{\text{inv}}$ spectrum
          & {25, 35, ..., 65} \\
$hmul_{[\text{K}^+, \text{K}^-\!]}^{\text{min}}$
          & Minimum multiplicity of hits in the CDC for a kaon candidate
          & \mbox{[27, 27]}, \newline \mbox{[30, 31]} \\
$p_{\text{K}^-}^{\text{max}}$ [GeV/$c$]
          & Maximum momentum of K$^-$ candidate & 0.60, 0.65 \\

  \end{tabular}
 \end{ruledtabular}
\end{table}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

\noindent Disc fragmentation through gravitational instability has been considered a potential formation channel for giant planets and very low mass stars for some time \citep{Cameron1978,Boss_science}.  Self-gravitating protostellar discs that are prone to  fragmentation produce objects with masses greater than a few Jupiter masses \citep{Forgan2011a} at semimajor axes larger than $\sim$ 30 AU \citep{Rafikov_05, Matzner_Levin_05, Whit_Stam_06,Mejia_3,Stamatellos2008, intro_hybrid, Clarke_09,Kratter2009,Vorobyov2010,Forgan2011a, Forgan2013, Tsukamoto2014}.  Once the conditions for fragmentation are met, the process occurs rapidly, typically on dynamical timescales.  Simulations of collapsing molecular clouds show that when discs fragment, it is typically within one free-fall time of the star's own birth \citep{Walch2009,Walch2010,collapses}, which is somewhat necessary given that self-gravitating discs produce strong, non-local angular momentum transport that can quickly drain the disc's mass and curtail fragmentation \citep{Harsono2011, Forgan2011}.

While much of the early work in disc instability theory focused on the fragmentation process, recent efforts have attempted to establish the final products of disc fragmentation - for example, the recently-coined tidal downsizing hypothesis \citep{Nayakshin2010b,Nayakshin2010,Nayakshin2010a} which builds on earlier work on fragment evolution \citep{Helled2008,Helled2008a,Boley2010b,Boley2011}.  Each fragment (hereafter embryo) contains a sample of the disc's population of dust grains.  It can add to this sample by accreting from the disc \citep{Helled2006} and grain size evolution can continue.  Grains will grow and sediment to the pressure maximum at the centre of the embryo, and can potentially form a core if the solids' density is sufficiently high.

During this process, the embryo's position in the disc is also changing. Migration due to disc torques can move the embryo inward, and can result in tidal stripping if the embryo's contraction rate is slow.  In extreme cases, this process can completely disrupt the embryo, bequeathing its evolved grain populations back to the disc \citep{Nayakshin2012} and potentially causing outburst events as the embryo's remains are accreted by the star \citep{Dunham2012,Nayakshin2012a, Nayakshin2013}.

Embryos can escape some of these dangers by undergoing a final collapse to a bound object, either due to the dissociation of hydrogen molecules as in the case of young protostars (cf \citealt{Masu_98}), or due to hydrostatic instability in the embryo's gaseous, metal-rich envelope around its newly forming core \citep{Nayakshin2014a}.  This may allow them to survive even at relatively low semimajor axes.

Recent population synthesis models \citep{TD_synthesis, Galvagni2013}, taking some or most of these effects into account, have extended the modelling timescales of fragmentation from a few tens of thousands of years up to around 1 Myr.  However, these timescales are still too short compared with observations of potential disc fragmentation descendants.  For example, the planetary system HR 8799 is often regarded as a candidate for disc fragmentation \citep{Nero2009}, with at least four $\sim 10 M_{\rm Jup}$ planets orbiting at 15-70 au from the host star.  However, the minimum age of the system is approximately 30 Myr \citep{Baines2012}.  If disc fragmentation modelling typically ends at the epoch of disc dispersal ($\sim$ 1 Myr) then how can we be confident that the systems produced remain in the same configuration?

The low eccentricities observed in the HR 8799 system appears to be an example of a system stabilised by a mean motion resonance \citep{Godziewski2014} but it is possible that disc fragments can be excited to high eccentricities.  Systems with multiple orbiting objects are likely to experience dynamical instabilities if the orbital separations become sufficiently low \citep{Chambers1996,Papaloizou2001,Chatterjee2008}.  Such systems allow planet pairs to undergo close encounters, which can result in planet-planet scattering.  These events excite orbital eccentricities, after which tidal circularisation could be responsible for short period planets such as Hot Jupiters \citep{Weidenschilling1996,Fabrycky2007}, as well as the ``tightly packed'' multiple planetary systems revealed by the Kepler Space Telescope \citep{Raymond2009a}.  In some cases, the eccentricity is excited beyond unity, ejecting one or more objects from the system \citep{Rasio1996}.   Planet-planet scattering has been invoked for core accretion theory to explain the eccentricity distribution in current exoplanet data \citep{Adams2003,Ford_and_Rasio_08, Juric2008}.  Given that multiple disc fragments can form in a single system, it is not unreasonable to expect that fragment-fragment scattering to also be important, even while disc gas is still present \citep[cf][]{Moeckel2008}.

We should therefore expect ejected disc fragments to be progenitors of a variety of substellar objects.  In particular, we should expect to see signatures in the free-floating planet population \citep{Sumi2011}.  The origin of these objects remains the subject of debate.  Planet-planet scattering appears to be insufficient to explain the entire free floating population \citep{Veras2012}, but it presumably plays an important role in shaping it. \emph{In situ} formation of substellar objects in molecular clouds is possible for objects of tens of Jupiter masses \citep{Hennebelle2008,Hennebelle2009,Strigari2012,Andre2012, Palau2014}, and tidal stripping of substellar companions while inside the parent cluster is likely to play some role \citep{Parker2012,Craig2013}\footnote{We should also note that stellar mass loss in the post-main sequence epoch could produce free-floating planets \citep{Mustill2013,Veras2013}, but we shall focus on younger systems in this work.}.

Forming more massive objects at larger semimajor axis increases the likelihood of stripping by the host cluster \citep{Boley2012}, and scattering events can be more energetic with more massive objects (although potentially less frequent due to larger orbital periods).  This would suggest that disc fragmentation may be extremely efficient at forming free floating planetary mass objects and field brown dwarfs.

In fact, \emph{if observations characterise the free-floating planet population to a sufficient level of detail and accuracy, then this will place important constraints on the fraction of protostellar discs that fragment}, a hitherto unknown quantity that will be essential for disc fragmentation theory to have any predictive power.

The current mass distribution of free-floating objects does appear to have a sharp cutoff near the theoretical prediction of minimum initial fragment mass, viz. a few Jupiter masses \citep{Sumi2011}, but this is probably due to observational limitations.  We should also note that re-analysis of the sample gives reduced numbers of objects classified as unbound \citep{Quanz2012}.

In any case, we should make predictions for future microlensing missions such as WFIRST \citep{Spergel2013}, which will be able to place tighter constraints on the free floating planet population.  Equally, the bound, large semimajor axis population of planets and brown dwarf companions will also be probed by WFIRST, providing further constraints on the frequency of disc fragmentation in the Milky Way.

In this work we attempt to ascertain the dynamical fate of disc fragments after the protostellar disc has dispersed.  We use data from \citet{TD_synthesis}'s population synthesis models of disc fragmentation as input to two separate studies of $N$-Body dynamics.  In the first, we consider isolated planetary systems, and carry out fragment-fragment scattering experiments to ascertain the subsequent distribution of orbital elements and ejection rates.  In the second, we place planetary systems within clusters to investigate orbital evolution and tidal stripping.  In section \ref{sec:methods} we briefly describe the population synthesis model and the $N$-Body experiments; in section \ref{sec:results} we describe the results, in section \ref{sec:discussion} we discuss the implications for exoplanet and free floating planet observations, and in section \ref{sec:conclusions} we summarise the work.

\section{Methods} \label{sec:methods}

\subsection{The Initial Dataset}

\noindent We take our sample of planetary systems formed by tidal downsizing from the population synthesis models of \citet{TD_synthesis}.   We summarise this model briefly here.

We generate a large number of protostellar disc systems, which are evolved using a standard $\alpha$-viscosity \citep{Shakura_Sunyaev_73, Rice_and_Armitage_09, Clarke_09}, where $\alpha$ represents a turbulent viscosity produced by self-gravity.  The disc masses are selected to ensure that this pseudo-viscous approximation is not invalidated by non-local angular momentum transport \citep{Lodato2005,Forgan2011}.  The transport of angular momentum due to self-gravity results in accretion of disc material onto the central star. The discs are also subject to mass loss by X-Ray photoevaporation driving winds \citep{Owen2010,Owen2011}.  The disc-to-star mass ratio is fixed, but the central star mass and the X-Ray luminosity are independently varied, producing a variety of disc lifetimes consistent with observations \citep{Haisch2001}.

In each of these protostellar discs, we determine sites of disc fragmentation by demanding that the Jeans mass inside a spiral arm perturbation be decreasing rapidly.  By calculating both the Jeans mass and its rate of change, we determine where a disc will fragment, and also the initial mass of the fragment \citep{Forgan2011,Forgan2013}.  The fragments are initially on circular orbits, and are separated by 1.5 to 3 Hill Radii.  The discs possess anywhere between 1 and 5 fragments.

The fragments are then evolved according to the semi-analytic framework of Nayakshin \citep{Nayakshin2010a,Nayakshin2010b,Nayakshin2010}.  The initial stages of fragment collapse are akin to those of the first cores (FCs) of the star formation process \citep{Masu_98}.  Grains grow and sediment once they reach a critical size to become sensitive to the local gas drag.  This sedimentation is suppressed by turbulence and by high-velocity collisions shattering grains.  Once the density of grains in the centre of the fragment exceeds the local gas density, the grains become self-gravitating and rapidly establish a core.  As the fragment continues to collapse, the central temperature rises.  Once this temperature exceeds $\sim$ 1600K, core formation is halted due to grain vapourisation.  Note that we do not include effects such as the core assisted gas capture instability recently described by \citet{Nayakshin2014a}.

During the evolution of the fragment and its solids population, it also migrates in the disc according to standard Type I/Type II migration timescales.  If the fragment radius exceeds its Hill Radius, then the upper layers of the fragment will be tidally stripped, and the mass of the fragment decreases.  While all fragments are initially several Jupiter masses, and hence firmly in the Type II regime, tidal stripping can move them to a Type I migration regime, causing the fragment's inward radial velocity to rapidly increase, and potentially completely destroying the fragment as the Hill Radius continues to shrink.  Fragment destruction is a common outcome from the population synthesis model - typically around 40\% of all fragments suffer this fate.  Note that this may be an underestimate of the true destruction rate, as more sophisticated migration modelling in self-gravitating discs suggest that migration timescales are shorter than the standard ones used here \citep{Baruteau2011}.

By repeating this process many times, we generate a population of objects.  In \cite{TD_synthesis} we ran four different iterations of the model to create four separate populations, investigating the effects of modifying various free parameters.  In this work, we use outputs from a single run of the model, corresponding to an opacity powerlaw index of 1, and ``standard'' migration (i.e., the migration timescale is not arbitrarily increased or decreased).  The initial distribution of object mass and semi-major axis can be found in Figure \ref{fig:ics}.

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{m_vs_a_pkap1_final.eps}
\caption{The mass-semimajor axis distribution retrieved from the population synthesis model of \citet{TD_synthesis}.  This dataset includes brown dwarfs (objects with masses greater than 13 $M_{\rm Jup}$), and giant planets with and without cores. \label{fig:ics}}
\end{center}
\end{figure}

\noindent

\subsection{Fragment-Fragment Scattering \label{sec:method_scattering}}

\noindent The population synthesis model generates individual systems with multiple objects, providing a large number of systems to integrate.  We sample 10,000 systems from the above dataset, with a total of 23,325 orbiting objects in the sample.  Each system has at least two bodies in orbit of the host star: 2-body star-planet systems are not investigated.

The population synthesis model gives a star mass, and orbiting object masses and semi-major axes.  It assumes that all orbits are circular, and it gives no information as to the relative phases of the orbits.  Given our ignorance of these parameters, we assume that the eccentricity is initially zero, all bodies are coplanar, and we initially assign orbital mean anomalies uniformly between $[0,2\pi]$.  The $N$-Body problem is solved using a 4th order Hermite integrator, with an adaptive global timestep.  The integrator conserves energy to 1 part in $10^6$ over the course of each simulation.  Note that we do not consider tidal evolution of the planet's orbit due to the star \citep[e.g.][]{Jackson2008}, and hence planets with highly eccentric orbits do not experience any eccentricity damping.  Collisions are not modelled in this analysis.

Each system is evolved for $10^6$ years.  This is still somewhat smaller than the observed ages of objects which we would like to compare to disc fragment models.  We select this relatively low simulation time to reduce computational expense and evolve more systems - typically, systems that produce scattering events express this instability within a few tens of thousands of years, which is to be expected given the objects' relatively close orbital spacing \citep{Chambers1996,Chatterjee2008}.  After this phase of instability, systems generally settle to quasistable orbital configurations.  Selecting this relatively short integration time will restrict our predictive power to some degree: for example, some systems will undergo eccentricity pumping that may lead to ejection on timescales greater than $10^6$ years.  In this sense we give conservative limits on both the eccentricity distribution and the ejection rate due to scattering.

\subsection{Simulating a Clustered Environment}{\label{sec:method:cluster}}

Most star-forming regions are observed to be filamentary \citep[e.g.][]{Arzoumanian2011}, which naturally leads to a hierarchically self-similar, or substructured spatial distribution of stars \citep{Cartwright2004,Sanchez2009,Gouliermis2014}. This substructure appears to be erased over time in observed star-forming regions \citep{Schmeja2008,Sanchez2009}, which is a consequence of 2-body relaxation \citep{Parker2012d,Parker2014b}.

\citet{Parker2014b} showed that dense regions that form with substructure and sub-virial or virial velocities will form a bound cluster \citep[see also][]{Allison2010}, whereas those with globally supervirial velocities (i.e.\,\,unbound) will retain substructure and possibly form large-scale association-like complexes.

The key parameter for assessing the impact of the star-forming environment on planetary systems is the initial local density. \citet{Parker2012} show that even a supervirial (unbound) region which is dense will alter the orbits of planetary systems.

It is important to note that a significant fraction of the dynamical interactions which will affect the planetary orbits occur in the early stages, even before the cluster has formed from the substructured star-forming region. For this reason, the effects of the star-forming environment will occur \emph{at the same time as the fragment-fragment scattering}.

Currently, simulations that simultaneously evolve multi-planetary systems within a clustered environment are in their infancy due to technical limitations, although recent progress has been made \citep{Hao2013,Liu2013,Pacucci2013}.

In this paper, we perform direct N-body star cluster simulations, in which the effects of interactions in the cluster on planetary-mass objects are explicitly computed (assuming the planet forms a binary system with its host star).  We randomly sample plantetary masses and semi-major axes from the distribution in \citet{TD_synthesis} before fragment--fragment scattering takes place, and determine the fraction that a) become free-floating, b) have their eccentricities significantly altered and c) have their semimajor axes significantly altered in the clustered environment. We then also use the simulations of \citet{Parker2012} with 1\,M$_{\rm Jup}$ single planets placed at 1, 5, 30 and 100\,au, which enables us to determine the fraction of systems that will be affected by the clustered environment which first have had their orbits altered by the fragment--fragment scattering experiments.

Our $N$-body simulations of star-forming regions are set-up in the same way as those in \citet{Parker2012}. The regions each have 750 stars drawn from  \citet{Kroupa2002} initial mass function of the form:
\begin{equation}
 N(M)   \propto  \left\{ \begin{array}{ll}
 M^{-1.3} \hspace{0.4cm} m_1 < M/{\rm M_\odot} \leq m_2   \,, \\
 M^{-2.3} \hspace{0.4cm} m_2 < M/{\rm M_\odot} \leq m_3   \,,
\end{array} \right.
\end{equation}
and we choose $m_1$ = 0.1\,M$_\odot$, $m_2$ = 0.5\,M$_\odot$, and  $m_3$ =
50\,M$_\odot$.

If a star has a mass of $>$2.5\,M$_\odot$ we pair it in a binary system with another star drawn from a flat mass ratio distribution \citep{Metchev2009,Reggiani2011,Reggiani2013} and a semimajor axis and orbital eccentricity drawn from the period distribution of binaries in the field \citep{Raghavan2010}\footnote{The binary properties of massive stars may actually differ signficantly from the Solar-mass stars in \citet{Raghavan2010}. However, the presence (or not) of binaries in a cluster does not change the resulting fraction of planets that are affected by dynamics \citep{Parker2012}, so the exact details of their set-up do not matter.}.  This is to prevent massive stars that may evolve significantly over the 10\,Myr timescale of the simulations from being paired with a planet. In the first set of simulations, all stars with lower masses are paired with a planet with mass and semimajor axis from the distribution in \citet{TD_synthesis}\footnote{In reality, many stars $<$2.5\,M$_\odot$ are also in binary systems. However, in order to fully sample the distribution from \citet{TD_synthesis} we pair each star with a planet, rather than a stellar companion.}, and in the second set of simulations, the planets all have a mass of 1\,M$_{\rm Jup}$ and are placed at 1, 5, 30 and 100\,au from the host star. The planets have zero eccentricity initially.

We note that this set-up may in fact over produce systems with planets formed via GI. However, the purpose of this paper is to perform a simple numerical experiment to assess the fractions of these planets that may then become free-floating due to interactions in a clustered environment. The initial stellar binary distribution that we adopt may not be an accurate representation of that in star-forming regions (which are largely unconstrained, e.g. \citealt{King2012}), however it is likely that the initial population is field-like \citep{Bate2012, Parker2014}. The deficit of stellar binaries in our calculations may underestimate the fraction of planets that become free-floating, as an disruptive interaction with a (relatively) massive binary system would be more likely than an interaction with a single star or star-planet system due to the higher collisional cross-section of that binary, although in earlier work \citet{Parker2012} do not find a strong dependence on stellar binarity.

The star-star and star-planet systems are then placed randomly in a substructured fractal distribution, as described in \citet{Goodwin2004a,Allison2010,Parker2014b}. A fractal distribution is the most convenient way of creating hierarchical substructure, because the level of substructure is quantified by a single number, the fractal dimension, $D$. Star-forming regions may, or may not be truly fractal \citep{Elmegreen2001}, but can usually be described by a fractal dimension \citep{Cartwright2004}.

Our fractal star-forming regions have a moderate amount of substructure $D = 2.0$, and a radius of 1\,pc. This results in a median density of $\sim$5000\,M$_\odot$\,pc$^{-3}$ in our simulations, which is somewhat higher than most local star-forming regions \citep{Bressert2010}. However, the `average' or dominant star-forming region which produces planet-hosting field stars is still unknown, and our simulations are consistent with the initial conditions which may have formed nearby star clusters such as the ONC \citep{Allison2010,Allison2011,Parker2011c}. For this reason, our simulations of the evolution of star-forming region probably represent an upper limit to the maximum density that planetary systems experience during their formation. The virial ratio of our star-forming regions is set to $\alpha_{\rm vir} = 0.3$, so that they are sub-virial and collapse to form a bound cluster after 1\,Myr.

The star-forming region are evolved for 10\,Myr using the \texttt{kira} integrator in the \texttt{Starlab} environment \citep[e.g.][]{Zwart1999,Zwart2001}, which utilises a 4$^{th}$-order Hermite scheme. We check whether a star--planet system is bound by determining the total energy of the system and the proximity of nearest neighbours \citep{Parker2009}. If the system has negative total energy, and the particles are mutual nearest neighbours, we consider it to be bound.  We do not include stellar evolution in the simulations.

\section{Results}\label{sec:results}

\subsection{Fragment-Fragment Scattering}\label{sec:results:frag}

\subsubsection{Properties of Bound Objects}

An example of a configuration that remains stable can be seen in Figure \ref{fig:stable}.  This system consists of a 1.18 $M_{\rm \odot}$ star, and 2 orbiting objects with masses of 14.6 $M_{\rm Jup}$ and 30.1 $M_{\rm Jup}$, and semimajor axes of 8.8 AU and 36.5 AU respectively.  This ordering of fragment masses - increasing mass with increasing separation - is a natural consequence of the Jeans mass formalism for initial fragment mass (see e.g. Figure 2 of \citealt{Forgan2011a}).  In this case, the mass ordering is preserved throughout the simulation, and both objects assume stable, low-eccentricity orbits (Figure \ref{fig:stable}).  This is primarily due to the more rapid orbital migration of the inner object separating the fragments by nearly 28 AU before self-interactions are activated.

\begin{figure*}
\begin{center}$
\begin{array}{cc}
\includegraphics[scale=0.4]{a_stable.eps} &
\includegraphics[scale=0.4]{e_stable.eps} \\
\end{array}$
\caption{Semimajor axis (left) and eccentricity (right) as a function of time for a single system that maintains a stable configuration.  The central star mass is 1.18 $M_{\rm \odot}$, and the 2 orbiting objects are 14.6 $M_{\rm Jup}$ (Object1) and 30.1 $M_{\rm Jup}$ (Object2) respectively.  \label{fig:stable}}
\end{center}
\end{figure*}

Figure \ref{fig:scattering} shows a system with a 0.8 $M_{\rm \odot}$ star and 4 orbiting objects, with masses 4.7 $M_{\rm Jup}$, 2.8 $M_{\rm Jup}$, 18.5 $M_{\rm Jup}$ and 7 $M_{\rm Jup}$ and semimajor axes of 12.2, 26.9, 47.6 and 55.9 AU respectively.  The tidal downsizing process does not preserve the mass ordering of fragments, resulting in an unstable configuration.  After several orbit crossings in the first few thousand years, the object second-closest to the star is ejected, and the eccentricities of the remaining objects are boosted.  In this relatively extreme case, the periastron radius of object number 3 is reduced to of order a solar radius, and would almost certainly undergo Roche lobe overflow and be completely destroyed and accreted.  Object 1 also comes within a relatively close distance, and would likely have its eccentricity damped by tidal evolution (as well as undergoing significant radiative heating and evaporation).

\begin{figure*}
\begin{center}$
\begin{array}{cc}
\includegraphics[scale=0.4]{r_scattering.eps} &
\includegraphics[scale=0.4]{rperi_scattering.eps} \\
\end{array}$
\caption{Distance from the star (left) and periastron radius (right) as a function of time for a single system that undergoes strong fragment-fragment scattering.  The central star is 0.8 $M_{\rm \odot}$, with the 4 orbiting objects having masses 4.7 $M_{\rm Jup}$ (Object1), 2.8 $M_{\rm Jup}$ (Object2), 18.5 $M_{\rm Jup}$ (Object3) and 7 $M_{\rm Jup}$ (Object4). \label{fig:scattering}}
\end{center}
\end{figure*}

\begin{figure*}
\begin{center}$
\begin{array}{cc}
\includegraphics[scale=0.4]{mvsa.eps} &
\includegraphics[scale=0.4]{evsa.eps} \\
\end{array}$
\caption{Mass versus semimajor axis after $10^6$ years of $N$-Body integration (left), and eccentricity versus semi-major axis after $N$-Body integration (right).  The initial eccentricity of all objects is zero.   \label{fig:mvsa_evsa}}
\end{center}
\end{figure*}

How are these processes encoded in the statistics of these bound objects? Figure \ref{fig:mvsa_evsa} shows the masses and eccentricities of the bound objects against semi-major axis at the end of the integration.  While the original distribution of objects with semimajor axis is largely preserved, with the median before and after integration being 29.9 and 24.2 AU respectively (see also Figure \ref{fig:ahist}), it is also apparent that scattering produces a long tail in the semimajor axis distribution out to $\sim 1 pc$, which skews the mean from 31.6 AU before integration to 2338.5 AU afterwards.

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{ahist.eps}
\caption{Semimajor axis distribution for all fragments that remain bound to their host star before (dashed green line) and after (solid blue line) $10^6$ years of $N$-Body integration. \label{fig:ahist}}
\end{center}
\end{figure}

The eccentricity of all objects is initially zero, as the population synthesis model is not equipped to evolve this quantity.  After integration, the objects occupy the entire range of eccentricities permitting bound orbits ($0 \leq e<1$).  The two bands in the right panel of Figure \ref{fig:mvsa_evsa} are indicative of the typical spacing between two fragments in a 3-body system, the most common outcome of the population synthesis model.  The dashed green line indicates a periastron radius of 0.1 AU - objects which fall to the left of this line can be expected to undergo tidal interactions to reduce their eccentricity (although we do not model this effect).  The red line indicates a periastron radius of $1R_{\odot}$ - on some occasions, objects achieve periastra close to this value, and are likely to undergo Roche lobe overflow and be partially or completely disrupted.

\subsubsection{Properties of Ejected Objects}

\noindent Of the 23,325 objects participating in these integrations, 5,956 achieve eccentricities greater than unity, and are ejected from the system.  Figure \ref{fig:mhist} compares the mass functions of all objects in the simulation (dashed green) to the objects ejected from their host systems (solid blue).  The general shape of both functions is similar - they share two peaks at $\sim 3 M_{\rm Jup}$ and $\sim 50 M_{\rm Jup}$.  However, it is clear that lower mass objects are preferentially ejected, enhancing the first peak and reducing the relative number of objects with masses greater than $\sim 20 M_{\rm Jup}$.  The origin of the two peaks is linked to the properties of the disc and the spacing between disc fragments.  As we have used G stars exclusively in this analysis, we only sample a small part of the mass function for disc fragment ejecta.  Increasing the mass of the star increase the typical fragment mass, however the stellar IMF will reduce their contribution to the final function.

Figure \ref{fig:ejectedmvsv} shows the distribution of ejecta velocities, as well as the relationship between mass and velocity.  The mean velocity of the ejecta is 2.4 $km s^{-1}$, with a median of 1.85 $km s^{-1}$.  There appears to be no correlation between the ejecta mass and their velocity.  This is in accord with ensembles of 3D hydrodynamical simulations of disc fragmentation \citep{Rice2003a,Stamatellos2009}.

As we will see in the next section, these velocities are similar to those found in the cluster experiments, and as such we cannot diagnose the formation mechanism of a free floating body from velocity alone.

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{ejectedmhist.eps}
\caption{Mass distribution of all fragments (dashed green line) and all ejected objects (solid blue line) after $10^6$ years of $N$-Body integration. \label{fig:mhist}}
\end{center}
\end{figure}

\begin{figure*}
\begin{center}$
\begin{array}{cc}
\includegraphics[scale=0.4]{ejectedvhist.eps} &
\includegraphics[scale=0.4]{ejectedmvsv.eps} \\
\end{array}$
\caption{The velocity distribution of ejected objects (left), and the mass-velocity distribution of ejecta (right). \label{fig:ejectedmvsv}}
\end{center}
\end{figure*}

\subsection{The Effect of Cluster Environment}

We begin by placing the population from \citet{TD_synthesis} into our subvirial (collapsing) star-forming regions. In Fig.~\ref{ecc_sma_regions} we show the orbitial eccentricity versus semimajor axis after 10\,Myr of evolution. The open circles are primordial, or birth systems -- the planet is still orbiting its parent star. The red plus signs represent captured systems -- the planets first become free-floating and then are captured into an orbit around a low- or intermediate-mass star ($m < 2.5$M$_\odot$) which no longer has its planet, or a more massive star which has lost its stellar partner and has accquired a planet. As in \citet{Parker2012}, the captured planets typically have high eccentricity and high inclinations (often retrograde orbits) and their semimajor axes are often 100s -- 1000s\,au \citep[see also][]{Perets2012}.

\begin{figure}
\begin{center}
\rotatebox{270}{\includegraphics[scale=0.4]{Plot_Or_SP_C_F2p0_S1FF10_ecc_sma.ps}}
\end{center}
\caption[bf]{Orbital eccentricity versus semimajor axis for planets after 10\,Myr of dynamical evolution in star-forming regions. We have summed together
the results of 10 simulations with the same initial conditions. Circles indicate `primordial' planetary systems
and crosses show captured planets.}
\label{ecc_sma_regions}
\end{figure}

The fraction of planets that become free-floating is $f_{\rm FFLOP} = 0.13 \pm 0.03$, the fraction that have their eccentricity raised to $e > 0.1$ is $f_{\rm e, alt} = 0.12 \pm 0.02$ and the fraction that have their semimajor axis altered by $\pm$10\,\% is $f_{\rm a, alt} = 0.057 \pm 0.009$. These fractions are not too dissimilar to those reported in \citet{Parker2012} for planets originally on 30\,au orbits, and the median semimajor axis in the distribution from \citep{TD_synthesis} is of order 30\,au.

We show the velocity distribution of the free-floating planets after 10\,Myr in Fig.~\ref{vel_regions}.  Panel a) gives the 3D velocity magnitude for planets still `observationally bound' to the cluster (within two half-mass radii, $2r_{1/2}$ of the centre), and panel b) gives the same property for the objects that are unbound (outside of two half-mass radii of the centre) in panel (b). We show the median velocity of the planets by the dashed line, and the median velocity of both stars and planets in the region by the dot-dashed line.

The free-floating planets that are still bound to the cluster have very similar velocities to the median, whereas those that are ejected from the cluster have much lower velocities than stars that are also ejected from the cluster. Again, this is similar to the result reported in \citet{Parker2012} for planets originally on 30\,au orbits.

\begin{figure*}
  \begin{center}
\setlength{\subfigcapskip}{10pt}
\hspace*{-1.5cm}\subfigure[Free-floating planets inside $2\,r_{1/2}$]{\rotatebox{270}{\includegraphics[scale=0.4]{Plot_Or_SP_C_F2p0_S1FF10_vel_bound.ps}}}
\hspace*{0.3cm}
\subfigure[Free-floating planets outside $2\,r_{1/2}$]{\rotatebox{270}{\includegraphics[scale=0.4]{Plot_Or_SP_C_F2p0_S1FF10_vel_esc.ps}}}
\caption[bf]{The distribution of velocities of free-floating planets
  that are outside two half-mass radii from the centre of the cluster
  ($2\,r_{1/2}$) after 10\,Myr. The median 3-D velocity for planets is shown by the dashed
  line whereas the median 3-D velocity for stars AND planets shown by dot-dashed line.}
\label{vel_regions}
  \end{center}
\end{figure*}

We have seen that many of the planets from the population synthesis in \citet{TD_synthesis} are ejected onto wide orbits through dynamical interactions in star-forming regions that subsequently form a star cluster. As discussed in Section~\ref{sec:method:cluster}, fragment--fragment scattering acts to increase the eccentricity and semimajor axis of the planets at the same time as the dynamical evolution of the star-forming regions. For this reason, it is difficult to disentangle the effects of fragment--fragment scattering from interactions in the star-forming region.

The probability of dynamical interactions in a star-forming region removing a free-floating planet \emph{after} fragment--fragment scattering onto a wider orbit can be approximated from the result of simulations in which we place single planets at 1, 5, 30 and 100\,au. We show the fraction of free-floating planets (the black plus signs), and the fractions of planets with significantly altered semimajor axes (red crosses) and eccentricities (blue asterisks), as a function of semimajor axis is shown in Fig.~\ref{fiducial_regions}.

The fraction of free-floating planets as a function of initial semimajor axis, $a$, can be approximated by a linear relation (the dotted black line), but a better fit is a power law where

\begin{equation}
f_{\rm FFLOP} = 0.027a^{0.46},
\end{equation}

\noindent meaning that 50\% of planets would be free-floating if their semimajor axis were greater than 570\,au, and all planets with $a > 8350$\,au would become free-floating in these dense regions.  Similarly,

\begin{equation}
f_{\rm e, alt} = 0.0072a^{0.79},
\end{equation}

meaning 50\% of planets would have $e>0.1$ if their semimajor axis were greater than 214\,au and,

\begin{equation}
f_{\rm a, alt} = 0.0042a^{0.75},
\end{equation}

\noindent and as such 50\,per cent of planets would have their semimajor axis altered by $a \pm $10\,per cent if their semimajor axis were initially 586\,au.

We note that the above fractions pertain only to planets formed from tidal downsizing, and do not include planets which may form on much closer (and hence more tightly bound) orbits.

It is unlikely that self-gravitating discs could maintain outer radii of several hundreds of AU due to truncation by stellar encounters \citep{Rosotti2014}, but the scattering simulations indicate that fragments can be boosted to orbits where the probability of becoming free-floating is significant. The mean semimajor axis of fragments still bound to the host after 1 Myr of fragment-fragment scattering is significantly larger than the 50\% threshold values quoted above (albeit the median remains much lower).

\begin{figure}
\begin{center}
\rotatebox{270}{\includegraphics[scale=0.4]{fiducial_systems.ps}}
\end{center}
\caption[bf]{Fraction of systems after 10\,Myr of dynamical evolution in star-forming regions that are free floating (the black plus signs), that have their eccentricity altered to above 0.1 (the blue asterisks) and those that have their semimajor axis altered by more than 10\,per cent (the red crosses), for planets initially at 1, 5, 30 and 100\,au from the host star. We also show fits to the points via the coloured lines; solid lines are power law fits, and dotted lines are linear fits.}
\label{fiducial_regions}
\end{figure}

Finally, we note that the mass of the planet is unimportant when determining the fraction of systems that are liberated form their host star in a clustered environment. The reason for this is that the planetary mass is almost negligible with respect to the host star, and any impulsive encounter that removes the planet is usually orders of magnitude more energetic than the mass-dependent binding energy of the system \citep{Fregeau2006,Parker2013b}.

\section{Discussion} \label{sec:discussion}

\subsection{Limitations of the Analysis}

\subsubsection{Fragment-Fragment Scattering}

\noindent Perhaps the most important outstanding issue with the isolated fragment-fragment scattering simulations is that fragment interactions are switched on after the disc is dissipated.  The population synthesis model makes a simple approximation of a single fragment's migration through the disc, thanks to the torques produced by the spiral arms that the fragment excites.  Multiple fragments will excite multiple spiral features, and hence the fragments' migration will depend on the morphology and evolution of these features.  We therefore cannot fully capture evolution such as resonant migration while the disc is present, which could increase the probability of mean motion resonances in these isolated systems.

We assume that no further evolution of the fragment's internal structure occurs during the $N$-Body integrations.  The grains trapped within the envelopes of the fragments may continue to grow and evolve long after the disc has dissipated.  For those objects which have internal temperatures surpassing the grain sublimation temperature, this may be a reasonable approximation.  However, it is clear that fragments undergoing close encounters may experience tidal forces strong enough to disrupt equilibrium density profiles.  In the extreme case, fragments that collide and merge may produce more massive objects, with the subsequent shocks affecting the object's chemistry and thermodynamic balance.  We do not model collisions in this work, and so this behaviour remains uninvestigated.

Even discounting collisions, close encounters between two fragments may also affect their geochemistry.  The integrations do not include the effects of tides being raised on the surface of the fragments, which may damp the fragment's eccentricity, and dissipate heat in its interior.  Close encounters will produce significant tides, and will likely cause significant melting or sintering of any solid core inside a fragment, which will have consequences for observed terrestrial planets if the envelope is stripped (cf \citealt{Nayakshin2014}).

\subsubsection{Cluster Simulations}

Our simulations of cluster environments and fragment evolution have several limiting factors. Firstly, the stars (and planets) are assumed to form instantaneously, which is obviously unrealistic. However, hydrodynamical calculations of star-formation are still limited by computing power, meaning that a large statistical sample of simulations cannot be achieved. They also typically only follow the first $\sim$Myr of star formation \citep[e.g.][]{Bate2012} and so cannot be evolved past the point at which dynamical interactions become negligible.

Ideally, we would also like to follow the effects of multi--fragment scattering \emph{and} the effects of stellar fly-bys simultaneously. Efforts have been made to implement this in full cluster simulations \citep[e.g.][]{Hao2013,Pacucci2013,Liu2013}, but in the meantime we can learn much from single-fragment evolution and also use collisional cross sections to make analytical estimates for the effect of the cluster environment on planetary systems \citep[e.g.][and many more]{Laughlin1998,Bonnell2001,Adams2006,Parker2012,Craig2013}.

Interestingly, the fraction of extrasolar planetary systems that could have been affected by dynamical interactions in the natal star-forming environment is almost totally unconstrained. In the local solar neighbourhood, \citet{Bressert2010} used the local surface density around YSOs to argue that most nearby star-forming regions are low-density $<$100\,stars\,pc$^{-2}$ and protoplanetary discs and systems are likely to be unaffected by dynamical evolution. However, it is unclear whether the \citet{Bressert2010} sample is representative of most star-forming events that produced the Galactic field population and even if it is, \citet{Parker2012d} argue that 50\,per cent of local star-forming regions could be dense enough to affect planetary systems.

The simulated star-forming regions presented here typically have maximum median densities of $\sim 5000$\,stars\,pc$^{-2}$, and we suggest that our simulations are probably an upper limit on the amount of dynamical evolution that may have affected young planetary systems.

\subsection{Implications for Observations}

\noindent Despite a great deal of dynamical evolution, the majority of objects formed through disc instability exist at large semimajor axes.  Even if the objects form with low eccentricity, this can be pumped up to high eccentricity through fragment-fragment scattering and dynamical interactions with the parent cluster.  This implies that current exoplanet detection missions remain relatively insensitive to disc instability objects, with the exception of microlensing and direct imaging surveys \citep[eg][]{Bonavita2014}.  While planets captured onto high inclination orbits are not exclusively disc fragments - for example, the Lidov-Kozai mechanism can produce high inclination orbits for a variety of planet parameters \citep{Naoz2013} -  detecting giant planets and brown dwarfs with large semimajor axes, eccentricities and inclinations are supporting evidence for their early formation while inside their parent cluster, a key attribute of gravitational instability.

Equally, the observed ejection rate of objects due to scattering, and the efficacy of stripping by a cluster for semimajor axes above $\sim 100$ AU, confirms that a large fraction of disc instability objects will exist as free-floating planets and field brown dwarfs.  The velocity dispersion of these objects will be relatively low compared to that of the parent cluster (at least initially).  These targets are also amenable to microlensing surveys, as well as infrared surveys of substellar objects \citep[eg][]{Best2013,Pinfield2013,Beichman2014}.

\section{Conclusions }\label{sec:conclusions}

\noindent We have investigated the dynamical fate of self-gravitating disc fragments once the disc has dissipated.  By carrying out separate $N$-Body integrations of a) isolated systems undergoing fragment-fragment scattering and b) systems embedded in substructured star-forming regions (which collapse to form clusters) we show that the descendants of disc fragments fill a large area of parameter space associated with the orbital elements.  While most objects retain semimajor axes and eccentricities that are similar to those of their formation, a significant fraction experience strong orbit modification, in some cases being ejected from their parent systems.  Some are recaptured on highly inclined, eccentric orbits.  Others are scattered inward toward the parent star, and are likely to undergo Roche lobe overflow and be destroyed.

By combining these calculations with the results of our population synthesis models \citep{TD_synthesis}, we can see that disc fragments share a variety of fates.  Over a third are completely destroyed by tidal disruption (see also \citealt{Zhu2012}). Of the surviving fragments, at least a quarter are ejected from the system - if fragment-fragment scattering and cluster dynamics were considered in the same integration, it is likely this figure would be higher.  The remainder tend to occupy orbits with semimajor axes greater than 10 AU, with eccentricities spanning all possible values for bound orbits.  Any planets that are to subsequently form in these  systems via core accretion must do so in the shadow of this dynamical evolution, and under the perturbative forces generated by the descendants of disc fragments.

Characterising the mass and velocity distributions of the free-floating planet population, especially in the 1-10 $M_{\rm Jup}$ range where tidal downsizing is favoured and core accretion disfavoured, will provide important constraints on the frequency of fragmenting systems, to complement direct imaging surveys for massive objects still bound to their parent star (e.g. \citealt{Bonavita2014, Bowler2014, Desidera2014}).  This work presents the first statistically significant sample of disc fragment properties post-ejection for observers to help estimate this frequency.  Any frequency derived from observations will depend sensitively on both the ejection rate and fragment destruction rate.  Further work must refine estimates for both of these properties.

\section*{Acknowledgments}

DF and KR acknowledge support from STFC grant ST/J001422/1.  DF also acknowledges support from the ECOGAL ERC Advanced Grant Programme.  RJP acknowledges support from the Royal Astronomical Society in the form of a research fellowship.

\bibliographystyle{mn2e}