OGLE-2017-BLG-1038: A Possible Brown-dwarf Binary Revealed by Spitzer Microlensing Parallax

We report the analysis of microlensing event OGLE-2017-BLG-1038, observed by the Optical Gravitational Lensing Experiment, Korean Microlensing Telescope Network, and Spitzer telescopes. The event is caused by a giant source star in the Galactic Bulge passing over a large resonant binary lens caustic. The availability of space-based data allows the full set of physical parameters to be calculated. However, there exists an eightfold degeneracy in the parallax measurement. The four best solutions correspond to very-low-mass binaries near ($M_1 = 170^{+40}_{-50} M_J$ and $M_2 = 110^{+20}_{-30} M_J$), or well below ($M_1 = 22.5^{+0.7}_{-0.4} M_J$ and $M_2 = 13.3^{+0.4}_{-0.3} M_J$) the boundary between stars and brown dwarfs. A conventional analysis, with scaled uncertainties for Spitzer data, implies a very-low-mass brown dwarf binary lens at a distance of 2 kpc. Compensating for systematic Spitzer errors using a Gaussian process model suggests that a higher mass M-dwarf binary at 6 kpc is equally likely. A Bayesian comparison based on a galactic model favors the larger-mass solutions. We demonstrate how this degeneracy can be resolved within the next ten years through infrared adaptive-optics imaging with a 40 m class telescope.


INTRODUCTION
Microlensing is a phenomenon in which the path of light emitted from a distant star (the source) is bent by a curve in space-time, caused by a massive object (the lens). If the source is approximately behind the lens, as seen by an observer, it brightens as unresolved images of the source are formed about the Einstein ring that has angular radius where π rel = au(D −1 L − D −1 S ), M is the mass of the lens system, D L and D S are the distance to the lens and source, respectively, and κ = 4G/(c 2 au) ∼ 8.14 mas/M .
For transient alignments, where the closest angular separation of the source and lens is on the order of θ E or smaller, photometric microlensing events can be observed as increasing and decreasing apparent brightness of the combination of the source star and unresolved neighbors, including the lens. Because only the source light is magnified, the luminosity of the lens system does not directly contribute to the event detection rate. As a result, microlensing is uniquely sensitive to the detection of low-mass, dim lenses such as brown dwarfs (BD; for example, Gould Paczyński (1986).
A limitation of the microlensing method is that, for most microlensing events, the light-curve model leaves a degeneracy between the mass and distance of the lens. This degeneracy can in principle be resolved either by measuring two other parameters (the Einstein radius θ E , and the microlens parallax π E ) or by separately observing the lens and source some years after the event in high-resolution images. While θ E has been measured for most planetary and binary events published to date, π E has not.
For events with an extremely dim lens, proper-motion measurement via late time imaging is not feasible at typical lens distances, given current observing capabilities. Breaking the mass-distance degeneracy for very faint lens systems thus requires a measurement of the microlens parallax. The spatial separation between observers required to detect parallax at a single epoch depends on characteristics of the microlensing event, such as the distance to the lens system and duration of the event. Because of the large separation between Earth and the Spitzer Space Telescope (located more than 1 au distant from Earth), microlensing observations from Spitzer, in conjunction with those from Earth, provide a reliable means of measuring parallax. This uniquely wide separation is what motivated the Spitzer microlensing project (Yee et al. 2015).
Microlensing has been used to discover 34 BDs from beyond the local regime (Chung et al. 2019). So far, this extended population has demonstrated unusual dynamics, such as an unexpected number of counter-rotating BDs (Chung et al. 2019; Shvartzvald et al. 2019Shvartzvald et al. , 2017. It is unclear to what degree these extreme kinematics are representative of the population as a whole.
BDs are stellar-like objects that are not massive enough to maintain a sufficient core temperature for main-sequence hydrogen fusion. Though the more massive BDs are capable of lithium fusion, and most BDs are capable of deuterium fusion, these processes do not provide sufficient heat to stop BDs from gradually cooling as they radiate the heat generated during their formation. As a result, they are very faint and become fainter as they age. Deuterium fusion occurs in objects with masses of approximately > 13 M J . This is often adopted as a criterion to distinguish BDs from planets; objects below this mass are planets, be they bound to a stellar object or free floating. However this mass definition is sometimes in conflict with the formation definition: BDs form like stars and planets form in circumstellar disks.
All but five of the microlensing BDs have been detected as binary systems. The number of BDs detected in binaries makes up an artificially high proportion of the total number of detections because binary events have more easily detected finite-source effects and therefore are more likely to have their associated masses calculated. Some of these have member masses at about the deuterium fusion limit (Choi et al. 2013;Han et al. 2017; Albrow et al. 2018), supporting the arguments of Grether & Lineweaver (2006) and Chabrier et al. (2014) for a mass overlap between the gas-giant planet and BD regimes. Deuterium fusion has become an insufficient metric for classification between BD and gas-giant planets. These populations have distinct formation histories, which, though difficult to infer, provide a more meaningful way to separate them in the mass-overlap region.
The upper BD cutoff is defined by sustained hydrogen fusion. Studies evaluating the hydrogen burning limit are summarised in Table 5 of Dieterich et al. (2018), from which we deduce that the BD upper limit is in the range of (∼ 70 − 95 M J ). This variance has a large dependence on chemical composition (e.g., Chabrier & Baraffe 1997). Forbes & Loeb (2019) investigate the idea of over-massive BDs. These are theoretically formed through Roche lobe overflow. The result is that, with only the mass information to draw from, this cutoff is vague.
Little is known about the very low mass end of the stellar initial mass functions (IMF). The empirical IMFs of Kroupa (2001), Chabrier (2005), Thies & Kroupa (2007); Thies & Kroupa (2008), and Kroupa et al. (2013) show disparity with the theoretical IMFs deduced from analytical descriptions of pre-stellar-cloud core distributions ( . Empirical IMFs usually require assumptions about age and metalicity in order to determine the IMF from an observed luminosity function. Observationally, measuring a mass function across the entire stellar mass range is challenging because sampling the upper mass range requires massive star clusters, and sampling the lower mass range requires nearby clusters. With the closest massive clusters at distances of a few kiloparsecs, observing both ends of the mass function in one star cluster is not currently possible photometrically (Elmegreen 2009). Wegg et al. (2017) shows one way in which microlensing surveys can be used to probe the IMF of the inner Milky Way, although this method used an existing dynamical model to infer the masses from the timescales (t E ) of ∼ 4000 events and therefore is not purely empirical. The timescales considered were 2 days < t E < 200 days, which relates to the mass via Currently, photometric surveys are only capable of probing relatively bright and very local populations of BDs. For example, Rosell et al. (2019) quote a distance limit in their Dark Energy Survey catalog of "beyond 400 pc". This selection bias in observability provides a limited view of BDs, in distance, mass and age. Further detections of very-low-mass objects in binary systems, will help to clarify our understanding of the dynamical properties of BD populations and the low-mass end of the IMF, because such systems are likely to have been formed as part of the very-low-mass end of the IMF, not like planets in a circumstellar disk.
The following sections in this paper describe our analysis of microlensing event OGLE-2017-BLG-1038 and how we determined this event to be a BD binary. §2 describes the observations made of this event, and the data-reduction methods used. §3 outlines our analysis of the ground-based data and resulting conclusions about source star characteristics. §4 details our analysis of the space-based, Spitzer data and our final modeling results. The corresponding physical parameters for our most likely models are calculated in §5. In §6 we compare the relative probabilities of our best model solutions and then we discuss, in §7, how different assumptions of the galactic model, as well as selection effects, may influence these probabilities. The Korean Microlensing Telescope Network (KMT-Net; Kim et al. 2016) also discovered this event as KMT-2017-BLG-0363 and observed it in the V and I bands. OGLE-2017-BLG-1038 was observed in two overlapping KMTNet search fields (BLG03 and BLG43), from each of the three KMTNet telescopes: Cerro Tololo Inter-American Observatory (KMT-C), South African Astronomical Observatory (KMT-S), and Siding Springs Observatory (KMT-A). This resulted in a cadence of ∼ 15 minutes between successive observations. The KMTNet observations were also primarily made in the I band. However, occasional V -band observations were made to provide color information. Therefore, 12 sets of KMT-Net light curves were obtained for this event.
The end of the event was also observed by the Spitzer Space Telescope Infrared Array Camera (IRAC; Fazio et al. 2004) instrument at an approximately 1 day cadence. While both the KMTNet and OGLE observations were made as part of regular survey operations, the Spitzer observations were scheduled for this event specifically as part of a program to enable space-parallax measurements for microlensing events (Calchi Novati et al. 2015a; Yee et al. 2015). This event was selected for Spitzer observations on 2017 June 13 (HJD' = 7918.11) and met the objective criteria on 2017 June 19 (HJD' = 7923.95). Both of these selections took place before the binary nature of the event was recognized, i.e., when it was still believed to be a point lens. Members of the Spitzer Team first noticed that the event was anomalous on 2017 June 20 (HJD' 7925.04).
Kinematic measurements from the source star in this event, as well as surrounding field stars, were obtained from Gaia Early Data Release 3 (Gaia Collaboration et al. 2020. The ground-based data were reduced using difference imaging (Tomaney & Crotts 1996;Alard & Lupton 1998) procedures. The OGLE images were reduced with their custom difference image procedures (see Wozniak 2000). The KMTNet light curves were extracted from the images using pyDIA (Albrow 2017) software, and the Spitzer light curve was extracted by the methods detailed in Calchi Novati et al. (2015b).

GROUND-BASED ANALYSIS
The light curve of this event (see Figure 1) has a triplepeaked perturbation over a 5 day period (2017 June [22][23][24][25][26][27] with the three peaks showing smoothed curves, indicative of a resolved source crossing a caustic. Caustics are features of a multiple-lens system. Therefore, we began our modeling with a binary-lens model, which we ultimately found was sufficient to describe the light curves for this event. The binary-lens model is parameterized by (s, q, ρ, u 0 , α, t 0 , t E ), where s is the angular separation of the two lens masses in units of θ E , q is the mass ratio of the lens objects, ρ is the source angular radius in units of θ E , u 0 is the closest line-of-sight point of approach to the lens center of mass made by the source in its relative trajectory (again in units of θ E ), and t 0 is the time at which this happens (|u 0 | = u (t 0 ), where u(t i ) is the position of the source, projected onto the lens plane, at a given time, (t), α is the angle of the projected rectilinear source trajectory relative to an axis that passes through the lens masses, and t E is the Einstein radius crossing time (the time the source takes to travel an angular distance of θ E ). For simplification, the motions in these models were considered from the reference frame of the lens system. This meant that, for modeling purposes, the relative velocities of any of the bodies involved were attributed to the "source velocity".
Our analysis of the ground-based light curves began by performing a grid search over a fixed resolution on s, q, u 0 , and α, using point-source approximations away from the caustics, for their computational speed, and convolved magnification maps in high-magnification regions, where finite-source effects were significant. The other model parameters were fitted by χ 2 minimization with ρ values found by interpolating between grid points with discrete convolutions. These calculations The best 20 grid solution regions were further investigated using the Emcee sampler (Foreman-Mackey et al. 2013). For this process we used the more accurate Image Centered Inverse RAy Shooting (ICIRAS) (Bennett was applied to the source in these calculations. Two of the regions converged to the same, and significantly most likely solution, while the next most likely solution had a ∆χ 2 of ∼ 110 000, before renormalization. The geometry of this static, ground-based solution is shown in Figure 2, and the magnification curve, with groundbased data, is shown in Figure 1. The fitted model parameters are displayed in Table 1 as the Static model. The solution corresponds to a source passing over the edges of a large resonant caustic. We note that this solution corresponds to small negative blending for three of the data sources, though this is a normal occurrence for microlensing photometry in a very crowded bulge field (Park et al. 2004), especially for dim lenses. Table  1 shows F B /F S for the OGLE source, which is within 2 σ of being positive. The source fluxes for each data set, were found from a linear fit; where F S is the source-star flux, F B is the blended flux 2 , A i is the magnification at time t i , and F i is the observed total flux at time t i . This solution to the static model was used to renormalize the ground-based data uncertainties (see Yee et al. 2012), and the solution was then allowed to reconverge.

Lens Orbital Motion or Ground-Based Parallax?
Although the peaks of the light curve are well fitted by this static-lens, rectilinear-source model, there is a region between dates 7915-7922 where the model systematically underpredicts the data ( Figure 3). In Figure 4 we show the cumulative χ 2 as a function of time for each individual data set. All curves show significant jumps near 7915-7922, indicating that there is a real missing feature in our static model. Higher-order effects are required for the model to provide a good description of these data.
Common high-order effects in microlensing light curves are orbital parallax (motion of Earth during an event) and orbital motion of the binary-lens system. A known degeneracy exists between these. Suspecting the significance of one or both of these higher-order effects, we added them to the generative model, both collectively and separately. We approximated the orbital motion of the lens objects by allowing α and s to vary linearly with time, adding the model parameterṡ α andṡ. Modeling the parallax effect requires the introduction of two new parameters, (π E,N , π E,E ), which are components of the vector π E , where π E = π rel θE , and its direction is that of the lens-source relative proper motion. The introduction of measurable parallax breaks the reflected symmetry of the source trajectory about the lens axis; a trajectory above the lens axis is not equivalent to a trajectory below the lens axis (except in the limit that the source lies exactly on the ecliptic). We therefore modeled both positive and negative u 0 solutions in which parallax was considered. For those solutions with both parallax and lens orbital motion, we calculate β (the ratio of the projected kinetic to potential energy of the lens; An et al. 2002;Dong et al. 2009), where values less than unity indicate a lens system consistent with a bound orbit; In our investigations of the significance of these two higher-order effects (Table 1), we find that, alone, lens orbital motion describes the static model discrepancies better than parallax. Including both higher-order effects yields only a minor χ 2 reduction compared with the purely lens-orbital-motion model, and the lens-orbital-motion parameters change very little. (The low β values for these models show that the implied orbits are bound.) Conversely, the posteriors of the parallax model change drastically when lens orbital motion is added. We therefore conclude lens orbital motion is well constrained and sufficient to describe the deviation on the static model from 7915-7922. This model is illustrated by the dotted lines in Figure 3.

Source Color
Color-magnitude Diagrams (CMDs) were created for each KMTNet observation site and field with I and V data (KMTC-03; Figure 5, KMTC-43, KMTS-03, KMTS-43, KMTA-03, and KMTA-43). We use the normal KMT practice of adopting magnitude zero points of I ZP = 28 and V ZP = 28.65. The source-star fluxes, obtained from fitting the magnification model to each light curve, were used to find the source star's position on the corresponding CMDs. The source fluxes for the highest likelihood solution (ground based) are given in Table 2.
The red clump in each CMD was centroid fitted, and acted as a calibration for obtaining the intrinsic colors and magnitudes of the field. The galactic bulge red clump can be used to calibrate the CMD because its intrinsic color and magnitude are known to high precision.

INCLUSION OF SATELLITE DATA
Having a Spitzer light curve for this event meant that, despite there being very inconclusive orbital parallax signals in the ground-based data, parallax could still be measured (Refsdal 1966). In this section we describe our analysis of the space-based Spitzer data using typical error renormalization methods, discuss concerns over systematics errors in the data, and present an alternate approach to coping with such systematics. Figure 6 shows the raw Spitzer data and a corresponding magnification curve from estimating F S = 56.1 (as is suggested by the color comparisons of §3.2), F B = 0, and adopting the ground-based model. In this figure, we can see a clear, decreasing signal that has ∆F > 30 Spitzer flux units. The Spitzer data are inconsistent with very small parallax, as the shape of the magnification curve is not well represented by the static ground-based model, and no alternative values of F S and F B could bring them into agreement. At the time of the first Spitzer observation, the ground-based light curve is still exiting the cusp while the Spitzer data are clearly not. This is strong evidence for a parallax effect. At the same time, the required magnification change as seen from Spitzer (∆A ∼ 1.6) indicates that the parallax cannot be too large.

Satellite Parallax Degeneracies
When viewed from Spitzer, the angular source trajectory across the lens plane is offset by a vector (∆β, ∆τ )/θ E , in directions (perpendicular, parallel) to D ⊥ , the separation between Spitzer and Earth pro-  Note-Those solutions indicated to by "LOM" refer to the models in which lens orbital motion was included. The source magnitude uses a zero point of IZP = 28. N is the total number of light-curve data points. Solutions with β < 1 are consistent with a bound orbit, but can only be calculated for models including both lens orbital motion and parallax. Note-These values were calculated using an orbiting, binary-lens model, for each of the ground-based sources. The Spitzer source flux is an estimate based on comparative CMDs between the Spitzer field and the KMTC-03 field.
jected onto the lens plane. This vector is related to the parallax measurement, but can be more useful in understanding the parallax likelihood space in comparison with the caustic diagram representation of the event.
The two parameters (∆β, ∆τ ) can be mapped onto π E,E and π E,N , via π E = au D ⊥ (∆τ, ∆β) . The parallel offset is simply In the case of a single lens, the perpendicular offset suffers from a four-fold satellite parallax degeneracy, due to the exact circular symmetry of the magnification field about the lens (Refsdal 1966), as illustrated in Gould (1994). (The sign convention we adopt here is that a positive value of u 0 indicates that, during its projected trajectory, the source approaches the lens center of mass on its right hand side.) In general, this fourfold degeneracy usually reduces to twofold with the addition of a second lens body, as the resulting caustic features break the symmetry of the magnification field. However, for binary-lens events in which the trajectory runs approximately parallel to the lens axis (such as the current case), trajectories reflected about the lens axis result in similar magnification curves, in which case the four-fold degeneracy is retained (Zhu et al. 2015).
A grid-search approach was used to determine the most likely parallax-solution regions. With the inclusion of space-based data, the two parallax parameters (π E,N and π E,E ) were added to the model.
When performing the parallax grid search, the ground-based model parameters (including lens orbital motion) were fixed, and a maximum-likelihood search was performed for the Spitzer light curve over a large range of discrete π E,N and π E,E values.
This grid search indicated that there were four solution regions for the given ground-based model, with the two outer regions having much higher likelihoods (i.e, lower χ 2 ) than the two inner regions (Figure 7a). These four solutions regions represent the ±u 0,Spitzer degenerate trajectories relating to two distinct solution families. We refer to these families as close (c) and wide (w). The four solutions regions result from only −u 0,Earth and indicate that, including the +u 0,Earth trajectory, we have an eightfold degeneracy for this particular geometry.
Because the Spitzer data only cover the falling part of the light curve and cover no caustic feature, the light curve alone does not contribute very strong constraints on the parallax measurement. We have thus implemented in the modeling an additional χ 2 penalty term that weighted the fit toward a source-flux ratio (between KMT-C03 and Spitzer L) matching that inferred by the calculated (I − L) 0 source color, found in §3.2. This color-constraint (Shin et al. 2017) term was of the form The constraint changed the likelihood space of the parallax model. The four solutions-regions from the unconstrained grid, remained as features in the constrained grid. However, the close set of solutions have more comparable likelihoods to the wide set than in the unconstrained grid.
When comparing Figure 7a and Figure 7b, the reason for the four lobes in the likelihood space becomes apparent. For this event, ∆β approximately aligns with π E,N and ∆τ with π E,E . Simplistically, changing ∆τ moves the Spitzer-data nodes backward or forward in time along the Spitzer trajectory, whereas ∆β shifts the "parallel" space-based trajectory closer to or farther away from the ground-based trajectory. The lobes and connective contours in Figure 7a result from solutions for which the Spitzer data hug the leftmost cusps of the caustic of Figure 7b. Figure 8 shows a more restricted view of the groundbased trajectory for this set of solutions, and the caustics at key epochs in the light curve, which change over the course of the event due to the orbital motion of the lenses.
Within each wide or close set, the pairs are the previously predicted ±u 0,Spitzer degenerate solutions. A further four degenerate solutions are obtained by reflecting all trajectories in Figure 7b (ground and Spitzer) about the lens axis.
The eight degenerate solution regions were further investigated using emcee with both ground-based and Spitzer data, renormalized errors, and both parallax and lens orbital motion included in the model. All model parameters were left free to evolve for all instances. Model parameters for the resulting solutions are given in Table 3. They are all somewhat similar in likelihood with an overall range in ∆χ 2 ≤ 87. The best solution found was the c -/+ geometry. All close solutions were favored over the wide by a margin of χ 2 w −χ 2 c ≥ 8.96. The nonfavored close solutions have a range 12.48 < ∆χ 2 < 28.06.

Spitzer Systematic Errors
Before we can have faith in these Spitzer parallax measurements, we must first address concerns of systematics in the Spitzer light curve.   where the flux levels were F S < 5 and thus F S ∼ F B , in which case systematics on the order of 1 could be considered fractionally significant.
We now consider whether systematics in the Spitzer data are significant for this event. The Spitzer magnification curve has a bump between t = 7936 and t = 7941 (corresponding to ∆F 5 Spitzer flux units; see Figure 6) that is not produced by any of our best generative model solutions that incorporate satellite parallax (Section 4.1). This implies a systematic error and demonstrates the scale to which we can expect them in this specific Spitzer data set; a few flux units over timescales of around 5 days. This is a higher ∆F perturbation than is expected for Spitzer systematic on a smooth curve (typically ∆F 1 Spitzer flux units) The parallax terms in the model are sensitive to small contiguous perturbations in the data, especially for those data after t = 7955, where flux changes of a few units change the shape of the slope enough to result in different parallax measurements, which affect the resulting physical solutions.   For this event, we have a Spitzer source flux much larger than the expected blend flux, a light curve with clearly and significantly decreasing flux, and baseline observations. Therefore, we would not ordinarily expect systematics to play a major role in this case. However, this event is somewhat sensitive to systematics in the baseline and shows evidence of similar systematics elsewhere in the light curve. We are therefore cautious of the effects systematic error in the Spitzer data may have on our conclusions.

Modeling Spitzer Errors
In an attempt to properly consider the apparent systematic errors in the Spitzer data, we have included in our model an error-bar renormalization parameter and two Gaussian process (GP) parameters.
Gaussian processes were first introduced in microlensing event analysis by Li et al. (2019). In this paper they used a GP to model source variability, rather than systematics, as well as a traditional inflated-error-bar scaling method. The GP method achieved better results in their case, as evidenced by the residuals in their Figure 1. However, they adopt their inflated-error-bar scaling model due to multiple practical and theoretical concerns. The practical issues they raise are how to cope with different blending effects between observation sources and how to perform error re-scaling. The blending issue is not relevant in our case because we only apply a GP model on the Spitzer data set. The theoretical issues they raise are in regards to choice of GP kernel and the possibility of degeneracies between the microlensing and GP parameters, for which they saw no evidence in their posterior distributions. We also saw no evidence of degeneracies between microlensing and GP parameters in our posterior distributions. In regards to the choice of GP kernel we tested both the exponential (described below) and Matern 3/2 kernels and found no significant difference between the results. We did not test the kernel used in Li et al. (2019) as it is meant for modeling quasi-periodic variations.
The degenerate solutions of Section 4.1 have reduced χ 2 values that imply that the Spitzer flux uncertainties have been underestimated by factors of between 2 and 5 times before renormalization. Because these factors change for each solution we include a multiplicative Spitzer error renormalization as a free parameter and consequently the likelihood must change to include the penalty ln P S = −N ln S, where S is the Spitzer error scaling factor and N refers here to the number of Spitzer data points.
Simultaneously we included an exponential GP model to fit the systematic features in the Spitzer light curve using the Celerite package (Foreman-Mackey et al. 2017). This replaces the vector of data variances with a data covariance matrix, We use a GP kernel where τ nm = |t n − t m |, and a and c are the GP model parameters.
The GP likelihood is then where r is the vector of (data -model) residuals.
The results of this modeling are displayed in Table  4. We find that inclusion of these three new model parameters has little effect on the microlensing parameters of all eight solutions, although the spread of likelihood values between solutions does change. With the GP parameters included in the model, our best solution is no longer the c -/+ but the c +/-, although, by a very small margin. The light curve for this model is shown in Figure 9. The full family of close solutions are all similarly likely, −2∆ ln L < 2.3, where we consider −2∆ ln L as an effective ∆χ 2 .

Model Comparison
When we compare the likelihoods using the standard analysis approach (Table 3) and GP analysis (Table 4), both favour the close solutions. For the close solutions, the range of ∆χ 2 < 28 using the standard approach and ∆χ 2 eff < 3 using GP. 3 The physical properties M tot and D L of all four close solutions are in agreement between the standard and GP approaches to within 1.2 σ. Therefore, the physical interpretation is the same in both cases. This is not true of all the degenerate wide-family solutions. While the large-parallax solutions remain most disfavoured between approaches, with matching M tot and D L values (masses at or below the deuterium fusion limit, 2 kpc away), the small-parallax solutions tell a different story. Using the standard approach, all wide solutions are disfavoured by a ∆χ 2 > 37. However, using the GP approach the small-parallax solutions have ∆χ 2 eff values of 7 and 15, within the ∆χ 2 range of close solutions using the standard approach. The physical interpretation is also different for these two solutions. The physical properties M tot and D L differ between approaches by < 4σ.
Our interpretation is that the physical solutions are not equally sensitive to systematic errors in the Spitzer data. The posteriors of π E,E are wider using the GP approach than the standard approach, especially the wide solutions for which the extrapolated trajectories do not cross caustics. It appears that the parallax measurement (particularly π E,E ) is proportionally more affected for smaller parallax solutions, making them more sensitive to systematic errors, but that the affect this has on the close solutions is limited by the nearness to a caustic crossing, which has a dominating effect on the likelihood space. Whether these conclusion are true in general is an interesting thought for future work.
While inflating error bars may be the correct approach for accommodating noise in data that is approximately Gaussian, it is appropriate to use a correlated noise approach where there are obvious systematic trends. The apparent perturbations in our Spitzer data are not represented by any of our best model solutions and therefore show that the errors in this data set are clearly correlated on time scales of a few days. However, the importance of using a correlated noise approach varies for our different solutions families and we believe that the importance of such modelling in other Spitzer events would also be dependent on many event-specific-properties.
Whether or not we consider the expense of a GP approach necessary, in our case, depends on the ∆χ 2 eff ranges we are prepared to accept. If we accept solutions at the ∆χ 2 eff 90 level, all eight degenerate solutions are valid, whether or not a GP is included. However, at the ∆χ 2 eff 50 level, we would reject the w -/+ and w +/solutions using the standard approach. Using the GP approach, we would accept all of these solutions, with w -/-and w +/+ converging into significantly different physical lens compositions.

Angular Einstein Radius
There exist empirical relations for determining the angular size of a star from its intrinsic color and magnitude. According to Kervella & Fouqué (2008), the most appropriate of these relations for non-M-type giants are those found in Nordgren et al. (2002) and van Belle (1999). 4 We use the Nordgren et al. (2002) surface brightness relation, specifically for non-variable giant stars (their Equation 12), Using the empirical color-color relations of Bessell & Brett (1988) for giant stars we find the (V − K) 0 equivalent of the intrinsic source color, (V − I) 0 , that was calculated from the CMDs, (V − K) S,O = 2.57 ± 0.09. The solutions for the models including higher-order effects have effectively identical θ * = 7.6 ± 0.5 µas.
θ E was calculated using the fitted ρ value for each solution, where θ E = θ * /ρ. The light-curve data provided good coverage of the caustic crossing and therefore ρ ±(4, 5) was well constrained, and almost identical, in our models. The calculated value of θ E for all solutions is θ E = 0.29 ± 0.02 mas.
Knowing θ E gives an angular scale to the geometric models.

Mass, Distance and Separation
The intrinsic I-band magnitude of the source star was previously calculated by comparing its fitted Iband magnitude to the mean red-clump magnitude on a CMD. By assuming the intrinsic red-clump magnitude and that the source star is at the distance of the average red-clump star in the CMD field, we find D S = 7.85 ± 0.06 kpc.
With values for D S , θ E and π E , the degeneracy in Equation 1 is broken, and the mass and distances can be calculated for each solution. Given the fitted parameters π E,E and π E,N , θ E , and D S , the distance to the lens was found, using where π rel = θ E π E . Knowing the distance to the lens system and θ E in angular units, the lens geometry can be calculated in absolute terms. The masses for each of the lens components, their projected separations, and the distances to the lens system are given in Tables  3 and 4. All large-parallax solutions, and both smallparallax wide solutions, are consistent with BD binary lenses of varying masses. However, the small-parallax close solutions are consistent with M-dwarf binaries, where the mass of the smaller of the binary objects (m 2 = 110 +20 −30 M J ) is very near the BD upper cut-off (∼ 70 − 95 M J ) and therefore may or may not be large enough for hydrogen fusion, depending mostly on its chemical composition.

Proper Motion and Velocity
The relative lens-source heliocentric proper motion was determined via  The source star for this event was observed by Gaia (EDR3 4063557344313009920), and hence its heliocentric proper motion is precisely measured as µ S,hel (N, E) = (−5.7, −7.7) ± (0.2, 0.3) mas yr −1, 5 relative to quasars in the distant universe. The source is ∼ 1σ due west of the centroid (see Figure 10). This means that a bulge lens is more easily accommodated, provided that direction of µ rel is roughly east. Similarly, the µ rel direction most consistent with a disk lens is northeast, although this direction is also very plausible for a bulge lens.
The heliocentric lens proper motion is calculated via The unexpected outcome of our µ L calculations is that none of the eight degenerate solutions align well with the disk or bulge dispersions, as shown in Figure 10. However, this demonstrates a misleading aspect of proper motion comparisons in that closer objects have higher proper-motions given the same tangential velocity. The lens proper motion relates to the heliocentric lens velocity via where distance is expressed in kiloparsecs, µ L,hel is in miliarcseconds per year, and 4.74 is a conversion factor so that v L,hel is in kilometers per second. These physical parameters for each solution can also be found in Tables 3 and 4. From Figure 10 we can see that the source is a fairly kinematically typical bulge star, lying on the 1σ contour of the Gaia field bulge dispersion.
Comparisons of the lens velocities, from each of the eight degenerate solutions, with disk and bulge dispersions from Gaia EDR3 are shown in Figure 11. These empirical dispersions are used for demonstrative purposes only. All eight lens solutions have unusual velocities when compared to typical disk stars, with the w +/+ and both +/-lens solutions rotating about the galactic center more slowly than typical disk stars, the w -/-and c +/+ counterrotating, and the -/+ and c -/-solutions seemingly moving through the disk, with large b velocities. The solutions are all less exceptional when compared with bulge kinematics, although only the small-parallax solutions have distances that allow for the lens to be a bulge member according to current galactic density models (e.g., Han & Gould 2003). The velocities of the w -/-, c +/+, w +/-, and c +/-solutions also appear consistent with the retrograde microlensing group. 5 Here we have doubled the published errors, as recommended by Rybizki et al. (2021).

SOLUTION PROBABILITIES
The somewhat uncommon physical parameters compel us to look at our solution probabilities more cautiously and holistically than a purely likelihood-based comparison. One problem with the likelihood calculation is that, formally, it relies on the assumption that our data are Gaussian distributed, with accurate uncertainties. Practically, this is never true for microlensing photometry. However, for this analysis, we apply Bayes theorem as though they were Gaussian.
The probability of a system having the solutionspecific proper motion or velocity, mass, and distance is also an important factor. We therefore calculate the probability factor ln z that determines the relative detection probability of each solution given a galactic model, with a bias to incorporate their relative light-curve-fit likelihoods.
We compute the galactic probability (Equation 15 of Gould (2020)) using a modified version of the Galactic Bayesian code described in Herrera-Martín et al. There is a common wisdom in microlensing analysis that small-parallax events are more probable than their large-parallax degenerate counterparts. This is known as the Rich argument, as detailed in Calchi Novati et al. (2015a). For single-lens events and binary-lens events for which the lens axis and source trajectory are approximately parallel (as in this case), if the true parallax solution is the smaller parallax solution it will always generate a large-parallax degenerate counterpart. The reverse, however, is not always true. The ratio of these probabilities (Rich factor) is implicitly accounted for in our galactic models (Gould 2020).
At the low galactic latitude of our event, and especially given the calculated distances to the lens of the large-parallax solutions, one would expect lens bodies to be members of the galactic disk. However, at a distance of ∼ 6kpc (as in the c -/-and c +/+ cases), it is possible that the lens is a member of the bulge population. Our galactic modeling of c +/+ showed that it is on the order of 100 times more likely to be a member of the bulge than the disk, whereas for c -/+ this was more like 1400 times more likely to be a member of the disk than the bulge. Currently, our galactic model most highly disfavours the counter-rotating BD solutions, with disklike distances (c +/-and w +/-with −2∆ ln z, without a light-cure likelihood bias, of 24.45 and 34.12, respectively) .
It is worth noting that ln z is based on a galactic model and therefore implicitly favors solutions matching our expectation of kinematic, mass, and density dispersions. Even the kinematic dispersions displayed in Figures 10  and 11 are informed by mostly bright stars and may not be truly representative of the dispersions of much dimmer objects, of which we know very little. Some healthy skepticism needs to exist around the model's completeness, especially considering the high proportion of microlensing BDs with unusual proper motions.
To determine how representative these retrograde detections are of the BD population as a whole, we must must first have a good understanding of the innate selection biases in microlensing events, for or against these extreme proper motions. However, if we were to downweight the light-curve likelihood based on the knowledge that our errors are not Gaussian, we will generally favour the low parallax solutions.

DISCUSSION
In our analysis of event OGLE-2017-BLG-1038, we fit a binary lens model including higher order effects: lens orbital motion and parallax. We include space-based data from Spitzer and model systematic errors in these data. We have a resulting eightfold solution degeneracy in this event. These solutions have total lens masses ranging from 0.027 − 0.27M . We also included in our probability comparison a galactic probability for each lens configuration. After these processes we find that our most probable solutions are the c +/+ and c -/-, both with masses of m 1 170 M J and m 2 = 110 +20 −30 M J (0.16 and 0.11 M ), separated by 1.7 au, at a distance of 6.0 kpc. The companion masses for these solutions are near the upper limit for BDs (the hydrogen burning limit). The lens systems for the c +/+ and c -/-solutions have tangential velocities of v L,hel (l, b) = (−358, −126) km s −1 and v L,hel (l, b) = (9, 113) km s −1 , respectively.
The c -/-solution has a minutely higher galactic probability than c +/+ with −2∆ ln L = 1.09. They are equally likely when considered in the context of both the light-curve fit and the galactic model.
Favouring these solutions over the large-parallax, close-family solutions (m 1 22.5 and m 2 13.7; D L = 2.33; v L,hel (l, b) = (−11, 88) km s −1 and v L,hel (l, b) = (−174, −21) km s −1 for c -/+ and c +/-, respectively) relies on our being confident in the galactic model for very-low-mass objects. Evidence from other microlensing events suggest that we do not understand the kinematic structure of BDs at distances of D < 4 kpc. To date, three BD systems have been discovered using microlensing that appear to be counterrotating with respect to the disk (Chung et al. 2019;Shvartzvald et al. 2019Shvartzvald et al. , 2017. These microlensing members lie very much in the plane of the disk and explanations for their characteristics, which we consider here, are that they are members of the disk with extreme motions; they are halo members with a coincidental disk alignment; they are members of a counterrotating population of verylow-mass objects (as suggested by Shvartzvald et al. 2019); or, they are evidence of an oversimplified galactic model. The physical parameters of the lens of this event raise the question as to whether or not OGLE-2017-BLG-1038 is another member of this group.
One explanation for extreme kinematics for a lowmass disk lens is that the disk may have a larger velocity dispersion for lower mass objects. If we assume that the lens was born in a cluster, it may have received a kick from an interaction with a star, and a binary will have a higher scattering cross section for such an interaction. Cluster dissolution has been extensively modeled  (2020) show most stars escaping with velocities < 10 km s −1 , relative to their parent cluster. It is therefore very unlikely that such an escapee would be travelling ∼ 100 km s −1 , or more, opposed to the disk. For globular clusters, higher mass objects preferentially wind up in tight binaries, whose members can be expelled at very high velocities (Hut et al. 1992a,b), but such expulsions are likely to account for a tiny fraction of all stars. This appears an unlikely origin for these counterrotating low-mass objects.
Another aspect of the galactic model that may be misunderstood is the bulge density model. We propose that a mass dependent spatial cutoff could explain the observed abundance of counterrotating BDs. If we consider that the bulge extends further for lower mass objects, then at D < 4 kpc the mass independent model would significantly underrepresent lower mass objects belonging to the bulge population and therefore having extreme (when compared to neighbouring disk stars) kinematics. Density models are fit to observational data and therefore are specifically fit to objects much larger than our inferred lens and those of the aforementioned retrograde BDs.
Another explanation may be that the lens is a halo star. Halo stars are known to have a much larger velocity dispersion, and their mean galactic rotation is much smaller than the disk (Du et al. 2018;Posti et al. 2018). While this large velocity dispersion could explain the kinematics of the other retrograde BD stars, it is a leap to make that assumption here, when it is not unlikely that the lens belongs to the bulge.
Are these retrograde BD detections the first members of a new class of object? At this stage, the characterization of these events as an independent population is speculative. Their existence as a discrete population affects the way we view the galactic probability of this solution, because such a population is not represented in the galactic model. Even if a misunderstood selection effect or aspect of the galactic model is responsible for their overabundance in detection, such an effect is not included in our current probability calculations. More needs to be known about this retrograde group before the significance of this solution can be truly understood.
The analysis of more low-mass lens events will provide new insights into the very-low-mass end of the mass function and its density and kinematics. There is little observational evidence to constrain any of these distributions at present. It is always possible that low-mass BDs are far more numerous than currently known and are currently represented by our galactic model. Whatever the case, for low-mass lenses, we believe that selection of a solution based on typical disk kinematic arguments is unlikely to be valid. The same reasoning leads us to believe that we cannot categorically claim this lens as a member of either a bulge, halo, or retrograde BD population. A more complex consideration of selection biases and possible population dynamics (beyond the scope of this paper) would be required.
A more empirical means of confirming the small parallax configuration would be to observe the lens photometrically. The hydrogen burning host and likely hydrogen-burning companion, corresponding to the small-parallax, close-family solutions, are bright enough to be visible at their implied lens distances (D L = 6 kpc). Given the relative proper motions of these solutions (µ rel,hel = 9.0 mas yr −1 ), we could expect the separation of source and lens to be sufficient for them to be resolved with the advent of infrared adaptive optics imaging from the coming generation of 40 m class telescopes. This is not true of the solutions near the planet-BD boundary, which are too dim to be resolved, no matter the angular separation between source and lens.
We expect first light for Multi-AO Imaging Camera for Deep Observations (MICADO) on the 39 m European Extremely Large Telescope (EELT) to be 2030. Kim et al. (2021) have argued, by scaling the work of Bowler et al. (2015) with the Keck coronograph, that an EELT coronograph could achieve ∆K = 11 contrast at 77 mas. By 2030 the angular separation of the lens and source will be ∼ 115 mas. Using the mass-luminosity function of Just et al. (2015) and the previously calculated source-star K magnitude, we estimate ∆K = 9.2 between the source star and the primary lens body for the M-dwarf solutions (c -/-and c +/+). Therefore the composition of this lens, be it BD or M-dwarf, can be verified with astrometric follow-up at the expected first light of MICADO on EELT.

SUMMARY
In this paper we report our analysis of microlensing event OGLE-2017-BLG-1038, with data from KMTNet, OGLE, and Spitzer. Ground-based data show the event is due to a giant source passing across a fold and cusp of a resonant caustic, due to a rotating binary lens. The analysis of the combined Spitzer, KMT, and OGLE light-curve data resulted in eight degenerate satelliteparallax solutions. With a GP model fit to the Spitzer data to account for systematic effects, the best solutions are the four belonging to the close family. Of these solutions the small-parallax solutions both have masses of M 1 170 +40 −50 M J (an M-dwarf) and m 2 = 110 +20 −30 M J (at the BD/M-dwarf cutoff). The large-parallax solutions are both comprised of a BD binary with m 1 = 22±2 M J and m 2 = 14 ± 1 M J . Inclusion of a detection probability based on a galactic model favors the small-parallax solutions. However, this approach to appraising solutions may be biased by an incomplete description of the distribution of very-low-mass objects in the galaxy and should not rule out solutions with similar light-curve-fit likelihoods. Late-time imaging could be used to reject these low-mass BD solutions, since an M dwarf should be visible given sufficient lens-source separation, but a low-mass BD binary will not.