High lithium-7 ($\mathrm{^7Li}$) abundances in giants are indicative of non-standard physical processes affecting the star. Mechanisms that could produce this signature include contamination from an external source, such as planets, or internal production and subsequent mixing to the stellar surface. However, distinguishing between different families of solutions has proven challenging, and there is no current consensus model that explains all the data. The lithium-6 ($\mathrm{^6Li}$) abundance may be a potentially important discriminant, as the relative $\mathrm{^6Li}$ and $\mathrm{^7Li}$ abundances are expected to be different if the enrichment were to come from internal production or from engulfment. In this work, we model the $\mathrm{^6Li}$ and $\mathrm{^7Li}$ abundances of different giants after the engulfment of a substellar mass companion. Given that $\mathrm{^6Li}$ is more strongly affected by Galactic chemical evolution than $\mathrm{^7Li}$, $\mathrm{^6Li}$ is not a good discriminant at low metallicities, where it is expected to be low in both star and planet. For modeled metallicities ([Fe/H]$>-0.5$), we use a "best case" initial $\mathrm{^6Li/^7Li}$ ratio equal to the solar value. $\mathrm{^6Li}$ increases significantly after the engulfment of a companion. However, at metallicities close to solar and higher, the $\mathrm{^6Li}$ signal does not last long in the stellar surface. As such, detection of surface $\mathrm{^6Li}$ in metal-rich red giants would most likely indicate the action of a mechanism for $\mathrm{^6Li}$-enrichment other than planet engulfment. At the same time, $\mathrm{^6Li}$ should not be used to reject the hypothesis of engulfment in a $\mathrm{^7Li}$-enriched giant or to support a particular $\mathrm{^7Li}$-enhancement mechanism.

The presence of Jupiter is crucial to the architecture of the Solar System and models underline this to be a generic feature of planetary systems. We find the detection of the difference between the position and motion recorded by the contemporary astrometric satellite Gaia and its precursor Hipparcos can be used to discover Jupiter-like planets. We illustrate how observations of the nearby star $\varepsilon$ Indi A giving astrometric and radial velocity data can be used to independently find the orbit of its suspected companion. The radial velocity and astrometric data provide complementary detections which allow for a much stronger solution than either technique would provide individually. We quantify $\varepsilon$ Indi A b as the closest Jupiter-like exoplanet with a mass of 3 $M_{Jup}$ on a slightly eccentric orbit with an orbital period of 45 yr. While other long-period exoplanets have been discovered, $\varepsilon$ Indi A b provides a well constrained mass and along with the well-studied brown dwarf binary in orbit around $\varepsilon$ Indi A means that the system provides a benchmark case for our understanding of the formation of gas giant planets and brown dwarfs.

In this work we use structural properties of Milky Way's outer halo ($R_G > 25\,\mathrm{kpc}$) satellites (dwarf spheroidal galaxies, ultra-faint dwarf galaxies and globular clusters) derived from deep, wide-field and homogeneous data, to present evidence of a correlation in the S\'ersic index v/s effective radius plane followed by a large fraction of outer halo globular clusters and satellite dwarf galaxies. We show that this correlation can be entirely reproduced by fitting empirical relations in the central surface brightness v/s absolute magnitude and S\'ersic index v/s absolute magnitude parameter spaces, and by assuming the existence of two types of outer halo globular clusters: one of high surface brightness (HSB group), with properties similar to inner halo clusters; and another of low surface brightness (LSB group), which share characteristics with dwarf spheroidal and ultra-faint dwarf galaxies. Given the similarities of LSB clusters with dwarf spheroidal and ultra-faint dwarf galaxies, we discuss the possibility that outer halo clusters also originated inside dark matter halos and that tidal forces from different galaxy host's potentials are responsible for the different properties between HSB and LSB clusters.

Stars of a common origin are thought to have similar, if not nearly identical, chemistry. Chemical tagging seeks to exploit this fact to identify Milky Way subpopulations through their unique chemical fingerprints. In this work, we compare the chemical abundances of dwarf stars in wide binaries to test the abundance consistency of stars of a common origin. Our sample of 31 wide binaries is identified from a catalog produced by cross-matching APOGEE stars with UCAC5 astrometry, and we confirm the fidelity of this sample with precision parallaxes from Gaia DR2. For as many as 14 separate elements, we compare the abundances between components of our wide binaries, finding they have very similar chemistry (typically within 0.1 dex). This level of consistency is more similar than can be expected from stars with different origins (which show typical abundance differences of 0.3-0.4 dex within our sample). For the best measured elements, Fe, Si, K, Ca, Mn, and Ni, these differences are reduced to 0.05-0.08 dex when selecting pairs of dwarf stars with similar temperatures. Our results suggest that APOGEE dwarf stars may currently be used for chemical tagging at the level of $\sim$0.1 dex or at the level of $\sim$0.05 dex when restricting for the best-measured elements in stars of similar temperatures. Larger wide binary catalogs may provide calibration sets, in complement to open cluster samples, for on-going spectroscopic surveys.

The LMC hosts a rich variety of star clusters seen in close projected proximity. Ages have been derived for few of them showing differences up to few million years, hinting at being binary star clusters. However, final confirmation needs to be done through spectroscopic analysis. Here we focus on the LMC cluster pair NGC2006-SL538 and aim to determine whether the star cluster pair is a bound entity (binary star cluster) or a chance alignment. Using the MIKE echelle spectrograph at LCO we have acquired integrated-light spectra for each cluster. We have measured radial velocities by two methods: a) direct line profile measurement yields v$_r=300.3\pm5\pm6$ km/s for NGC2006 and $v_r=310.2\pm4\pm6$ km/s for SL538. b) By comparing observed spectra with synthetic bootstrapped spectra yielding $v_r=311.0\pm0.6$ km/s for NGC2006 and $v_r=309.4\pm0.5$ km/s for SL538. Finally when spectra are directly compared, we find a ${\Delta}v=1.08\pm0.47$ km/s. Full-spectrum SED fits reveal that the stellar population ages lie in the range 13-21 Myr with a metallicity of Z=0.008. We find indications for differences in the chemical abundance patterns as revealed by the helium absorption lines between the two clusters. The dynamical analysis shows that the two clusters are likely to merge within the next $\sim$150 Myr. The NGC2006-SL538 cluster pair shows radial velocities, stellar population and dynamical parameters consistent with a gravitational bound entity. We conclude that this is a genuine binary cluster pair, and we propose that their differences in ages and stellar population chemistry is most likely due to variances in their chemical enrichment history within their environment. We suggest that their formation may have taken place in a loosely bound star-formation complex which saw initial fragmentation but then had its clusters become a gravitationally bound pair by tidal capture.

Samples of reliably identified halo wide binaries are scarce. If reasonably free from selection effects and with a small degree of contamination by chance alignments, these wide binaries become a powerful dynamical tool, having provided one of the very few experiments capable of constraining the nature of dark matter in the Milky Way halo. Currently, however, the best available sample of halo wide binaries is confined to the solar neighborhood, and is plagued by small number statistics at the widest separations. We present the results of a program aimed to probe the wide binary population of the Galactic halo at significantly larger distances, and which informs future searches that could improve the statistics by orders of magnitude. Halo stars were taken from the Sloan Digital Sky Survey after analysing the Galactic orbits of stars in the reduced proper motion diagram. We then select candidate binaries by searching for pairs with small differences in proper motion and small projected separation on the sky. Using medium-resolution spectroscopy, a subsample of candidates is validated via radial velocities, finding that about 68\% of candidate pairs up to 20$^{\prime\prime}$ separation can be considered genuine halo wide binaries, to the limits of the available data. Precise parallaxes from \textit{Gaia} confirm that most of our selected pairs have their components at the same distances, independently confirming the robustness of our selection method. These results should prove valuable to guide the optimal assembly of larger catalogs of halo wide binaries from upcoming large databases, e.g., \textit{Gaia} and LSST.

We perform a search for wide binaries in the Kepler field with the prospect of providing new constraints for gyrochronology. First, we construct our base catalog by compiling astrometry for the stars observed by Kepler, and supplement it with parallaxes, radial velocities (RVs), and metallicities. We then mine our base catalog for wide binary candidates by matching the stars' proper motions, as well as parallaxes, RVs, and metallicities, if available. We mitigate the presence of chance alignments among our candidates by performing a comprehensive data-based contamination analysis in the proper motion versus angular separation phase space. Our final sample contains 55 binary candidates. A crossmatch of our pairs with the Second Data Release (DR2) from Gaia validates our candidates and confirms the reliability of our search method, particularly for $\varpi \gtrsim 2$ mas. Due to the implicit Kepler selection function and image scale per pixel, our binary search is incomplete for angular separations of less than 20 arcsec. We crossmatch our candidates with rotation period and asteroseismic ages catalogs, and find that our binary candidates do not follow a simple period-color relation, in agreement with previous studies. Two pairs have an age estimate for one component star and rotation period for its companion, positioning them as potentially new gyrochronology constraints at old ages. This is the first study that uses RVs and metallicities as criteria, rather than as a confirmation, in a binary search.

Featuring greatly improved astrometry and precise radial velocities (RVs), the second data release (DR2) from Gaia affords us a unique opportunity to test the validity of samples of stellar wide binaries. In a previous work, we presented a set of wide binaries identified in the joint Tycho-Gaia Astrometric Solution catalog. Here, we use DR2 to confirm the low contamination rate of this sample, thereby verifying the effectiveness of our algorithm for application to DR2 data.

We address the existence and origin of the lithium (Li) desert, a region in the Li - T_eff plane sparsely populated by stars. Here we analyze some of the explanations that have been suggested for this region, including mixing in the late main sequence, a Li dip origin for stars with low Li abundances in the region, and a possible relation with the presence of planets. To study the Li desert, we measured the atmospheric parameters and Li abundance of 227 late-F dwarfs and subgiants, chosen to be in the T_eff range of the desert and without previous Li abundance measurements. Subsequently, we complemented those with literature data to obtain a homogeneous catalog of 2318 stars, for which we compute masses and ages. We characterize stars surrounding the region of the Li desert. We conclude that stars with low Li abundances below the desert are more massive and more evolved than stars above the desert. Given the unexpected presence of low Li abundance stars in this effective temperature range, we concentrate on finding their origin. We conclude that these stars with low Li abundance do not evolve from stars above the desert: at a given mass, stars with low Li (i.e., below the desert) are more metal-poor. Instead, we suggest that stars below the Li desert are consistent with having evolved from the Li dip, discarding the need to invoke additional mixing to explain this feature. Thus, stars below the Li desert are not peculiar and are only distinguished from other subgiants evolved from the Li dip in that their combination of atmospheric parameters locates them in a range of effective temperatures where otherwise only high Li abundance stars would be found (i.e., stars above the desert).

From our recent catalog based on the first Gaia data release (TGAS), we select wide binaries in which both stars have been observed by the Radial Velocity Experiment (RAVE) or the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST). Using RAVE and LAMOST metallicities and RAVE Mg, Al, Si, Ti, and Fe abundances, we find that the differences in the metallicities and elemental abundances of components of wide binaries are consistent with being due to observational uncertainties, in agreement with previous results for smaller and more restricted samples. The metallicity and elemental abundance consistency between wide binary components presented in this work confirms their common origin and bolsters the status of wide binaries as "mini-open clusters." Furthermore, this is evidence that wide binaries are effectively co-eval and co-chemical, supporting their use for e.g., constraining age-activity-rotation relations, the initial-final mass relation for white dwarfs, and M-dwarf metallicity indicators. Additionally, we demonstrate that the common proper motion, common parallax pairs in TGAS with the most extreme separations (s $\gtrsim$ 0.1 pc) typically have inconsistent metallicities, radial velocities or both and are therefore likely to be predominantly comprised of random alignments of unassociated stars with similar astrometry, in agreement with our previous results. Finally, we propose that wide binaries form an ideal data set with which to test chemical tagging as a method to identify stars of common origin, particularly because the stars in wide binaries span a wide range of metallicities, much wider than that spanned by nearby open clusters.

We mine the Tycho-{\it Gaia} astrometric solution (TGAS) catalog for wide stellar binaries by matching positions, proper motions, and astrometric parallaxes. We separate genuine binaries from unassociated stellar pairs through a Bayesian formulation that includes correlated uncertainties in the proper motions and parallaxes. Rather than relying on assumptions about the structure of the Galaxy, we calculate Bayesian priors and likelihoods based on the nature of Keplerian orbits and the TGAS catalog itself. We calibrate our method using radial velocity measurements and obtain 6196 high-confidence candidate wide binaries with projected separations $s\lesssim1$ pc. The normalization of this distribution suggests that at least 0.6\% of TGAS stars have an associated, distant TGAS companion in a wide binary. We demonstrate that {\it Gaia}'s astrometry is precise enough that it can detect projected orbital velocities in wide binaries with orbital periods as large as 10$^6$ yr. For pairs with $s\ \lesssim\ 4\times10^4$~AU, characterization of random alignments indicate our contamination to be $\approx$5\%. For $s \lesssim 5\times10^3$~AU, our distribution is consistent with \"{O}pik's Law. At larger separations, the distribution is steeper and consistent with a power-law $P(s)\propto s^{-1.6}$; there is no evidence in our data of any bimodality in this distribution for $s \lesssim$ 1 pc. Using radial velocities, we demonstrate that at large separations, i.e., of order $s \sim$ 1 pc and beyond, any potential sample of genuine wide binaries in TGAS cannot be easily distinguished from ionized former wide binaries, moving groups, or contamination from randomly aligned stars.

The Gaia-ESO survey recently reported on a large sample of lithium (Li) abundance determinations for evolved stars in the rich open cluster Trumpler 20. They argue for a scenario where virtually all stars experience post main sequence mixing and Li is preserved in only two objects. We present an alternate explanation, where Li is normal in the vast majority of cluster stars and anomalously high in these two cases. We demonstrate that the Li upper limits in the red giants can be explained with a combination of main sequence depletion and standard dredge-up, and that they are close to the detected levels in other systems of similar age. In our framework, the two detected giants are anomalously Li-rich, and we propose that both could have been produced by the engulfment of a substellar mass companion of 16(+6-10) M_J. This would imply that ~5% of 1.8 solar mass stars in this system, and by extension elsewhere, should have substellar mass companions of high mass that could be engulfed at some point in their lifetimes. We discuss future tests that could confirm or refute this scenario.

A small fraction of red giants are known to be lithium (Li) rich, in contradiction with expectations from stellar evolutionary theory. A possible explanation for these atypical giants is the engulfment of a Li-rich planet or brown dwarf by the star. In this work, we model the evolution of Li abundance in canonical red giants including the accretion of a sub-stellar mass companion. We consider a wide range of stellar and companion masses, Li abundances, stellar metallicities, and planetary orbital periods. Based on our calculations, companions with masses lower than 15 M_J dissolve in the convective envelope and can induce Li enrichment in regimes where extra mixing does not operate. Our models indicate that the accretion of a substellar companion can explain abundances up to A(Li)~2.2, setting an upper limit for Li-rich giants formed by this mechanism. Giants with higher abundances need another mechanism to be explained. For reasonable planetary distributions, we predict the Li abundance distribution of low-mass giants undergoing planet engulfment, finding that between 1% to 3% of them should have A(Li)>1.5. We show that depending on the stellar mass range, this traditional definition of Li-rich giants is misleading, as isolated massive stars would be considered anomalous while giants engulfing a companion would be set aside, flagged as normal. We explore the detectability of companion engulfment, finding that planets with masses higher than ~7 M_J produce a distinct signature, and that descendants of stars originating in the Li-dip and low luminosity red giants are ideal tests of this channel.

We report the discovery of BD+20594b, a high density sub-Neptune exoplanet, made using photometry from Campaign 4 of the two-wheeled Kepler (K2) mission, ground-based radial velocity follow-up from HARPS and high resolution lucky and adaptive optics imaging obtained using AstraLux and MagAO, respectively. The host star is a bright ($V=11.04$, $K_s = 9.37$), slightly metal poor ([Fe/H]$=-0.15\pm 0.05$ dex) solar analogue located at $152.1^{+9.7}_{-7.4}$ pc from Earth, for which we find a radius of $R_*=0.928^{+0.055}_{-0.040}R_\odot$ and a mass of $M_* = 0.961^{+0.032}_{-0.029}M_\odot$. A joint analysis of the K2 photometry and HARPS radial velocities reveal that the planet is in a $\approx 42$ day orbit around its host star, has a radius of $2.23^{+0.14}_{-0.11}R_\oplus$, and a mass of $16.3^{+6.0}_{-6.1}M_\oplus$. Although the data at hand puts the planet in the region of the mass-radius diagram where we could expect planets with a pure rock (i.e. magnesium silicate) composition using two-layer models (i.e., between rock/iron and rock/ice compositions), we discuss more realistic three-layer composition models which can explain the high density of the discovered exoplanet. The fact that the planet lies in the boundary between "possibly rocky" and "non-rocky" exoplanets, makes it an interesting planet for future RV follow-up.

We carry out controlled $N$-body simulations that follow the dynamical evolution of binary stars in the dark matter (DM) haloes of ultra-faint dwarf spheroidals (dSphs). We find that wide binaries with semi-major axes $a\gtrsim a_t$ tend to be quickly disrupted by the tidal field of the halo. In smooth potentials the truncation scale, $a_t$, is mainly governed by (i) the mass enclosed within the dwarf half-light radius ($R_h$) and (ii) the slope of the DM halo profile at $R\approx R_h$, and is largely independent of the initial eccentricity distribution of the binary systems and the anisotropy of the stellar orbits about the galactic potential. For the reported velocity dispersion and half-light radius of Segue I, the closest ultra-faint, our models predict $a_t$ values that are a factor 2--3 smaller in cuspy haloes than in any of the cored models considered here. Using mock observations of Segue I we show that measuring the projected two-point correlation function of stellar pairs with sub-arcsecond resolution may provide a useful tool to constrain the amount and distribution of DM in the smallest and most DM-dominated galaxies.

Gas-rich galaxies in dense environments such as galaxy clusters and massive groups are affected by a number of possible types of interactions with the cluster environment, which make their evolution radically different than that of field galaxies. The dIrr galaxy NGC 1427A, presently infalling towards the core of the Fornax galaxy cluster, offers a unique opportunity to study those processes in a level of detail not possible to achieve for galaxies at higher redshifts. Using HST/ACS and auxiliary VLT/FORS ground-based observations, we study the properties of the most recent episodes of star formation in this gas-rich galaxy, the only one of its type near the core of the Fornax cluster. We study the structural and photometric properties of young star cluster complexes in NGC 1427A, identifying 12 bright such complexes with exceptionally blue colors. The comparison of our broadband near-UV/optical photometry with simple stellar population models yields ages below ~4x10^6 yr and stellar masses from a few thousand up to ~3x10^4 Msun, slightly dependent on the assumption of cluster metallicity and IMF. Their grouping is consistent with hierarchical and fractal star cluster formation. We use deep Ha imaging data to determine the current Star Formation Rate (SFR) in NGC 1427A and estimate the ratio, Gamma, of star formation occurring in these star cluster complexes to that in the entire galaxy. We find Gamma to be among the largest such values available in the literature, consistent with starburst galaxies. Thus a large fraction of the current star formation in NGC 1427A is occurring in star clusters, with the peculiar spatial arrangement of such complexes strongly hinting at the possibility that the starburst is being triggered by the passage of the galaxy through the cluster environment.

The Kapteyn moving group has been postulated as tidal debris from $\omega$ Centauri. If true, members of the group should show some of the chemical abundance patterns known for stars in the cluster. We present an optical and near-infrared high-resolution, high-S/N spectroscopic study of 14 stars of the Kapteyn group, plus 10 additional stars (the $\omega$ Cen-group) that, while not listed as members of the Kapteyn group as originally defined, have been nevertheless associated dynamically with $\omega$ Centauri. Abundances for Na, O, Mg, Al, Ca and Ba were derived from the optical spectra, while the strength of the chromospheric He I 10830 {\AA} line is studied as a possible helium abundance indicator. The resulting Na-O and Mg-Al patterns for stars of the combined Kapteyn and $\omega$ Cen-group samples do not resemble those of $\omega$ Centauri, and are not different from those of field stars of the Galactic halo. The distribution of equivalent widths of the He I 10830 {\AA} line is consistent with that found among non-active field stars. Therefore, no evidence is found for second-generation stars within our samples, which most likely rules out a globular-cluster origin. Moreover, no hint of the unique Ba-overabundance at the metal-rich end, well-established for $\omega$ Centauri stars, is seen among stars of the combined samples. Because this specific Ba pattern is present in $\omega$ Centauri irrespective of stellar generation, this would rule out the possibility that our entire sample might be composed of only first generation stars from the cluster. Finally, for the stars of the Kapteyn group, the possibility of an origin in the hypothetical $\omega$ Centauri's parent galaxy is disfavored by the different run of $\alpha$-elements with metallicity between our targets and stars from present-day dwarf galaxies.

We present the first study of high-precision internal proper motions (PMs) in a large sample of globular clusters, based on Hubble Space Telescope (HST) data obtained over the past decade with the ACS/WFC, ACS/HRC, and WFC3/UVIS instruments. We determine PMs for over 1.3 million stars in the central regions of 22 clusters, with a median number of ~60,000 stars per cluster. These PMs have the potential to significantly advance our understanding of the internal kinematics of globular clusters by extending past line-of-sight (LOS) velocity measurements to two- or three-dimensional velocities, lower stellar masses, and larger sample sizes. We describe the reduction pipeline that we developed to derive homogeneous PMs from the very heterogeneous archival data. We demonstrate the quality of the measurements through extensive Monte-Carlo simulations. We also discuss the PM errors introduced by various systematic effects, and the techniques that we have developed to correct or remove them to the extent possible. We provide in electronic form the catalog for NGC 7078 (M 15), which consists of 77,837 stars in the central 2.4 arcmin. We validate the catalog by comparison with existing PM measurements and LOS velocities, and use it to study the dependence of the velocity dispersion on radius, stellar magnitude (or mass) along the main sequence, and direction in the plane of the sky (radial/tangential). Subsequent papers in this series will explore a range of applications in globular-cluster science, and will also present the PM catalogs for the other sample clusters.

We present a program designed to obtain age-rotation measurements of solar-type dwarfs to be used in the calibration of gyrochronology relations at ages of several Gyr. This is a region of parameter space crucial for the large-scale study of the Milky Way, and where the only constraint available today is that provided by the Sun. Our program takes advantage of a set of wide binaries selected so that one component is an evolved star and the other is a main-sequence star of FGK type. In this way, we obtain the age of the system from the evolved star, while the rotational properties of the main sequence component provides the information relevant for gyrochronology regarding the spin-down of solar-type stars. By mining currently available catalogs of wide binaries, we assemble a sample of 37 pairs well positioned for our purposes: 19 with turnoff or subgiant primaries, and 18 with white dwarf components. Using high-resolution optical spectroscopy, we measure precise stellar parameters for a subset of 15 of the pairs with turnoff/subgiant components, and use these to derive isochronal ages for the corresponding systems. Ages for 16 of the 18 pairs with white dwarf components are taken from the literature. The ages of this initial sample of 31 wide binaries range from 1 to 9 Gyr, with precisions better than 20% for almost half of these systems. When combined with measurements of the rotation period of their main sequence components, these wide binary systems would potentially provide a similar number of points useful for the calibration of gyrochronology relations at very old ages.

Schwarzschild's orbit-superposition technique is the most developed and well-tested method available for constraining the detailed mass distributions of equilibrium stellar systems. Here I provide a very short overview of the method and its existing implementations, and briefly discuss their viability as a tool for modeling the Galaxy using Gaia data.

Establishing or ruling out, either through solid mass measurements or upper limits, the presence of intermediate-mass black holes (IMBHs) at the centers of star clusters would profoundly impact our understanding of problems ranging from the formation and long-term dynamical evolution of stellar systems, to the nature of the seeds and the growth mechanisms of supermassive black holes. While there are sound theoretical arguments both for and against their presence in today's clusters, observational studies have so far not yielded truly conclusive IMBH detections nor upper limits. We argue that the most promising approach to solving this issue is provided by the combination of measurements of the proper motions of stars at the centers of Galactic globular clusters and dynamical models able to take full advantage of this type of data set. We present a program based on HST observations and recently developed tools for dynamical analysis designed to do just that.

(ABRIDGED) We present a new Schwarzschild orbit-superposition code designed to model discrete datasets composed of velocities of individual kinematic tracers in a dynamical system. This constitutes an extension of previous implementations that can only address continuous data in the form of (the moments of) velocity distributions, thus avoiding potentially important losses of information due to data binning. Furthermore, the code can handle any combination of available velocity components, i.e., only line-of-sight velocities, only proper motions, or a combination of both. It can also handle a combination of discrete and continuous data. The code finds the distribution function (DF, a function of the three integrals of motion E, Lz, and I3) that best reproduces the available kinematic and photometric observations in a given axisymmetric gravitational potential. The fully numerical approach ensures considerable freedom on the form of the DF f(E,Lz,I3). This allows a very general modeling of the orbital structure, thus avoiding restrictive assumptions about the degree of (an)isotropy of the orbits. We describe the implementation of the discrete code and present a series of tests of its performance based on the modeling of simulated datasets generated from a known DF. We find that the discrete Schwarzschild code recovers the original orbital structure, M/L ratios, and inclination of the input datasets to satisfactory accuracy, as quantified by various statistics. The code will be valuable, e.g., for modeling stellar motions in Galactic globular clusters, and those of individual stars, planetary nebulae, or globular clusters in nearby galaxies. This can shed new light on the total mass distributions of these systems, with central black holes and dark matter halos being of particular interest.