#  Time-dependent potential impurity in topological insulator

Saurabh Pradhan  Department of Physics and Astronomy, Box 516, 75120, Uppsala
University, Uppsala, Sweden  Institut für Experimentalphysik, Freie
Universität Berlin, Arnimallee 14, 14195 Berlin, Germany  Jonas Fransson
Department of Physics and Astronomy, Box 516, 75120, Uppsala University,
Uppsala, Sweden

###### Abstract

We consider periodically driven potential impurities coupled to the surface
states of a two-dimensional topological insulator. The problem is addressed by
means of two models, out which the first model is an effective continuum
Hamiltonian for the surface states, whereas the Kane-Mele lattice model is our
second approach. While both models result in drastic changes in the local
density of electron states with increasing amplitude and frequency of the
driving field, the linearly low energy local density of electron states
remains in the continuum model, however, with an increased Fermi velocity. The
spectrum of the continuum model remains gapless under the emergence of new
impurity resonances near the Fermi energy. The Kane-Mele lattice model
represents a finite size system, with edge states appearing at the boundary of
the system. We, thus, consider the impurity at two different positions, one at
the boundary and one at the center of the lattice. In the former case, a
reduction and broadening of the low energy local density of electron states
result with increasing amplitude of the driving field. On the other hand,
there are no new resonances emerging in the spectrum. In the latter case, the
spectrum is gapped both in the absence of the impurity as well as for weak
amplitudes of the driving field, while the gap tends to fill up with impurity
states with increasing amplitude.

##  I Introduction

Classical and quantum systems show new phases with unexpected properties when
driven by an external time-dependent periodic field. In many cases, these
phases do not have any counterpart in static systems. For example, periodic
driving in a one dimensional systems generates chaotic motion of particles
chaos1 ; chaos2 ; chaos3 ; chaos4 . In molecular systems, temporal changes of
microscopic parameters or bias voltage may lead to various exotic non-
equilibrium states Fainshtein . The ratio between kinetic and potential energy
as well as lattice spacing can be varied in optical lattices to create non-
trivial states coldatom1 ; coldatom2 . Recently, periodically driven systems
have given rise new phases drive1 ; drive2 ; drive3 ; drive4 ; drive5 ; drive6
; drive7 ; drive8 ; drive9 ; drive10 , such as, Floquet block-states in
topological insulator ftimat1 ; ftimat2 , optical lattices ftiopt , and cold-
atom fticold  systems. Topological insulators have metallic edge states due to
the spin-orbit coupling and non-trivial topology of the band structure, while
sustaining an insulating gap in the bulk. Material such as Bi2Se3 and Bi2Si3
shows topologically protected Dirac cone, resulting from the properties of
these edge states ti1 ; ti2 ; ti3 . Moreover, due to spin-momentum locking,
back-scattering off potential impurities is prevented. In other words, the
definite chirality of the edge states does not allow transitions between the
states $|{\bf k}\uparrow\rangle$ and $|-{\bf k}\downarrow\rangle$ assisted by
scattering off time-reversal symmetric (non-magnetic) impurities.

The robustness of these edge mode has been proven theoretically by examining
the backscattering due to potential impurity in continuum and lattice model
edge1 ; edge2 ; edge3 ; edge4 ; edge5 ; edge6 . The protection of the Dirac
surface states against disorder is also seen in various experiments edge7 ;
edge8 ; edge9 . Magnetic impurities, on the other hand, do not keep these edge
modes intact which results in imperfect quantization of the conductance mag-
im1 ; mag-im2 ; mag-im3 .

Even though there are many theoretical and experimental studies of the effects
of static impurities on topological insulator, the effects of dynamical
impurities have not been considered previously. Dynamical local interaction
can give rise to various non-linear phenomena such as multi-photon
dissociation or excitation of atom or molecule when exposed under strong laser
field lesser1 ; lesser2 ; lesser3 . In this article, we consider an atomically
sharp time-dependent impurity potential, simulating the conditions of a local,
yet extremely focused, monochromatic laser field, affecting a single lattice
site. We show the emergence of impurity resonances, in accordance with
previously discussed for impurity resonances for stationary conditions edge4 ;
edge5 ; edge6 ; PhysRevB.94.075401 , which give rise to strong modifications
of the local density of electron states near the impurity site, with an
increasing number of resonance features with increasing amplitude and
frequency of the driving field. While the metallic edge states tend to develop
a sharp density peak around the Fermi level for impurities interacting with
the edge states, the density gap of the bulk states become increasingly filled
up with increasing scattering potential strength, for impurities directly
perturbing the bulk sites.

##  II Model and Floquet-Green Function

In this section, we describe the general formalism of the Floquet Green
function method. The Hamiltonian of the full quantum system is a periodic
function in time, $H(t)=H(t+\tau)$, where $\tau=2\pi/\Omega$ is the period of
the external driving field. As a direct consequence of the periodicity, we can
use the Floquet theorem Floquet1 ; Floquet2 ; Floquet3  which can be regarded
as a time domain equivalent to the Bloch theorem. Due to the explicit time-
dependence of the Hamiltonian, electrons can be excited to different energy
states. However, if the time-dependent potential is characterized by a single
frequency ($\Omega$), the energy difference between the final and initial
state should be an integer multiple of $\Omega$. This restriction gives rise
to an energy space representation of the Hamiltonian, the corresponding Green
function, as well as of operators.

We define the Floquet Hamiltonian, a Hermitian operator, for the generic time-
dependent Hamiltonian $H(x,t)$

$\displaystyle H^{F}(x,t)=H(x,t)-i\hbar\frac{\partial}{\partial t},$ (1)

which gives the Floquet green function written as

$\displaystyle\Bigl{(}\epsilon-H^{F}(x,t)\Bigr{)}G^{F}(x,x^{\prime};t,t^{\prime})=\delta(x-x^{\prime})\delta(t-t^{\prime}).$
(2)

The time periodicity of the Hamiltonian is inherited in the Floquet Green
function, which becomes periodic in both $t$ and $t^{\prime}$. Hence, we can
Fourier expand both the Green function and Hamiltonian into

$\displaystyle H^{F}(x,t)=$
$\displaystyle\sum_{\gamma}H^{F}_{\gamma}(x)e^{i\epsilon_{\gamma}t},$ (3a)
$\displaystyle G^{F}(x,x^{\prime};t,t^{\prime})=$
$\displaystyle\sum_{\alpha,\beta}G_{\alpha\beta}^{F}(x,x^{\prime})e^{i\epsilon_{\alpha}t-i\epsilon_{\beta}t^{\prime}}.$
(3b)

Here, the quasi-energy $\epsilon_{\alpha}$ is a conserved quantity. In this
way, a time-periodic driven system is reduced to an algebraic matrix equation.
Although the Fourier expansion reduces the complexity of the problem, the
dimensions of the corresponding Hilbert space is infinite. Therefore, there is
no exact analytical closed form of the Green function. In fact, even the exact
solution of a two-level system driven with linearly polarized light is not
known shirley .

Due to this intrinsic complexity, we approach the Floquet Green function
numerically, by considering a harmonic monochromatic time-dependent driving
field. The Hamiltonian, then, assumes the form

$\displaystyle\mathbf{H}=\mathbf{H}^{0}+2\mathbf{V}\cos\Omega t,$ (4)

where $\mathbf{H}^{0}$ is the time-independent part of the Hamiltonian whereas
$\mathbf{V}$ represents the coupling to the time-dependent driving field. The
corresponding Floquet Hamiltonian for the assumed driving field has block tri-
diagonal structure, according to

$\displaystyle
H^{F}_{\alpha\beta}=(\mathbf{H}^{0}-\alpha\hbar\Omega)\delta_{\alpha\beta}+\mathbf{V}(\delta_{\alpha+1\beta}+\delta_{\alpha-1\beta})$
(5)

The Floquet green function becomes

$\displaystyle(\mathbf{1}E_{\alpha}-\mathbf{H}^{0})G_{\alpha\beta}-\mathbf{V}(G_{\alpha+1\beta}+G_{\alpha-1\beta})=\mathbf{1}\delta_{\alpha\beta},$
(6)

where $E_{\alpha}=\epsilon+\alpha\hbar\Omega$. This matrix equation can be
solved iteratively by using the _matrix continued fraction_ method martinez1 ;
martinez2 . As a result, we obtain a recursive matrix equation for the Green
function

$\displaystyle\Big{[}\mathbf{1}E_{\alpha}-\mathbf{H}^{0}-\mathbf{V}_{\text{eff}}(E_{\alpha})\Big{]}G_{\alpha,\alpha}=$
$\displaystyle\mathbf{1},$ (7a)
$\displaystyle\mathbf{V}_{\text{eff}}(E_{\alpha})=$
$\displaystyle\mathbf{V}^{+}_{\text{eff}}(E_{\alpha})+\mathbf{V}^{-}_{\text{eff}}(E_{\alpha}).$
(7b)

Here, the effective potential is given by

$\displaystyle\mathbf{V}^{\pm}_{\text{eff}}(E_{\alpha})=$
$\displaystyle\mathbf{V}\frac{1}{\mathbf{C}_{\alpha\pm
1}-\mathbf{V}\frac{1}{\mathbf{C}_{\alpha\pm
2}-\mathbf{V}\frac{1}{\vdots}\mathbf{V}}\mathbf{V}}\mathbf{V},$ (8)

where $\mathbf{C}_{\alpha}=\mathbf{1}E_{\alpha}-\mathbf{H}^{0}$. Assuming that
the impurity sit at the origin $(\mathbf{r}=0)$, the effective potential can
be written

$\displaystyle\mathbf{V}^{\pm}_{\text{eff}}(E_{\alpha})=$
$\displaystyle|0\rangle\langle 0|\frac{1}{\mathbf{C}_{\alpha\pm
1}-|0\rangle\langle 0|\frac{1}{\mathbf{C}_{\alpha\pm
2}-\frac{1}{\vdots}}|0\rangle\langle 0|}|0\rangle\langle 0|$ $\displaystyle=$
$\displaystyle|0\rangle\langle 0|\frac{1}{\mathbf{C}_{\alpha\pm
1}-|0\rangle\langle 0|\mathbf{V}^{\pm}_{\text{eff}}(E_{\alpha\pm
1})|0\rangle\langle 0|}|0\rangle\langle 0|$ $\displaystyle=$ $\displaystyle
V^{\pm}_{\text{eff}}(E_{\alpha})|0\rangle\langle 0|.$ (9)

In general, we compute these effective potentials numerically by setting a
maximum frequency $E_{M}$ above which the effective potentials are zero, that
is, $\mathbf{V}^{\pm}_{\text{eff}}(E_{m})=0$, for all $E_{m}>E_{M}$. Typically
we take $m$ to be on the order of 100, which normally is sufficient for
convergence.

##  III Topological insulator surface

First, we will investigate the effect of time-dependent impurity on the edge
of three-dimensional topological insulators. The low energy effective model
within $k\cdot p$ approximation can be described by

$\displaystyle
H=v\sum_{k}\Psi_{\mathbf{k}}^{\dagger}[\mathbf{k}\times\vec{e}_{3}]\cdot\sigma\Psi_{\mathbf{k}}$
(10)

We are interested in the local properties of this system such as density of
states $N({\bf r},\omega)$, which is related to the (retarded) Green function
through the identity $N({\bf r},\omega)=-{\rm tr}\;{\rm Im}{\bf G}^{r}({\bf
r},{\bf r};\omega)/\pi$. As a function of the effective potential, we write
the Green function in momentum space as

$\displaystyle G_{k,k^{\prime}}=$
$\displaystyle\delta_{kk^{\prime}}G_{k}^{0}+G^{0}_{k}\frac{V_{\text{eff}}}{1-V_{\text{eff}}G^{0}_{0}}G^{0}_{k^{\prime}},$
(11)

where $G$, $G^{0}$, and $V$ are implicit functions of the energy.

We consider the driving potential
$\mathbf{V}=2A\delta(\mathbf{r}-\mathbf{r}_{0})\sigma^{0}\cos\Omega t$ with
frequency $\Omega$, amplitude $A$, at the position $\mathbf{r}_{0}$. In Fig.
1, we plot the local density of electron states for increasing impurity
amplitude $A$, panels (a) through (d). In absence of the impurity, Fig. 1 (a),
we retain the typical Dirac-like density of states with vanishing density of
state at the Fermi energy, as expected. However, the linear density of states
around the Fermi energy is preserved also at finite driving amplitudes, Fig. 1
(b) – (d). In addition to the linear low energy spectrum, high energy features
emerge with increasing $A$. These features, which appear symmetrically around
the Dirac point, are direct consequences of the excitations that are generated
by the time-dependent potential, much in analogy with the impurity side
resonances that are caused by vibrational defects on Dirac materials
PhysRevLett.110.026802 ; PhysRevB.87.245404 . In fact, the symmetric
appearance of the excitations is caused by the combination of positive and
negative scattering potentials ${\bf V}^{+}_{\text{efft}}$ and ${\bf
V}^{-}_{\text{efft}}$, respectively, each of which is responsible for the
equally distributed set of excitations on either the valence _or_ the
conduction side of the electronic band structure. Moreover, the number of
excitations increases substantially with increasing $A$, which we understand
to be an effect the increased order of $\mathbf{C}_{\alpha}$ that contributes
to the effective potential $\mathbf{V}_{\text{eff}}$. Therefore, the Green
function picks up an increasing number of higher energy modes. We also notice
that the bandwidth of the local density of states increases the stronger the
driving force. We refer this to the increasing number of contributing
excitations that become available for in the scattering processes and which
necessitates a redistribution of the total density onto an increased set of
excitations.

Figure 1: Local density of electron states at the impurity site with driving
frequency $\Omega=0.50$ for different values of the driving amplitude $A=$ 0.0
(a), 0.5 (b), 1.0 (c), 2.0 (d).

##  IV Kane-Mele model

In the second part of this article, we will consider a lattice, Kane-Mele,
model for a topological insulator. The model pertains to a tight-binding
Hamiltonian on a honeycomb lattice which is a straight forward generalization
of the Haldane model haldane . Here, we write

$\displaystyle H=$ $\displaystyle-t\sum_{\langle
ij\rangle\sigma}c^{\dagger}_{i\sigma}c_{j\sigma}+i\lambda_{so}\sum_{\langle\langle
ij\rangle\rangle\sigma}c^{\dagger}_{i\sigma}\sigma^{z}c_{j\sigma}+\mu\sum_{i\sigma}c^{\dagger}_{i\sigma}c_{i\sigma},$
(12)

where $t$ is the nearest-neighbor (NN) hopping integral, $\lambda$ is the
spin-orbit coupling which act as a next-nearest-neighbor hopping element,
whereas $\mu$ is the chemical potential which fixes the number of particle of
the system. Here, we consider a half-filled system for which $\mu=0$. Note
that spin-orbit coupling act as a purely imaginary hopping integral and differ
by a sign for up and down spin component. We set $t$ as our absolute energy
scale, and we have imposed open boundary conditions in the $x$-direction and
periodic boundary conditions along the $y$-direction. We calculate
$V_{\text{eff}}$ iteratively for a lattice size $100\times 50$. We consider
the same driving potential as above,
$\mathbf{V}=\delta(\mathbf{r}-\mathbf{r}_{0})2A\sigma^{0}\cos\Omega t$, with
the frequency $\Omega$, amplitude $A$ at $\mathbf{r}_{0}$.

Figure 2: (a) Total density of state of the Kane-Mele model for the
topological insulator with a single impurity at the middle of the lattice. The
vertical lines indicate the energy at which local density of state is plotted
in the bottom panel. Real space map of the local density of state for
different values of energy E=0.0 (b), 0.5 (c), 1.0 (d).

The plots in Fig. 2 display the total density of states of the Kane-Mele model
with a single impurity located at the center of the lattice. The red vertical
line signifies the energy values at which the spacial maps of the local
density of state are plotted in panels (b) through (d). The density of states
for this tight binding model on a honeycomb lattice vanishes at the Fermi
energy and increases linearly away from the Fermi energy. A finite spin-orbit
coupling opens up a gap in the spectrum. Due to the non-trivial topological
nature of this system, this gap only appears in the bulk while metallic edge
states appear at the boundary. In Fig. 2 (b), (c), this can be seen as an
enhanced density of states at the edges $x=0$ and $x=100$. Moreover, one
should notice the finite density of states around the impurity, a density
which oscillates and decays far away from it. At higher energies, however,
bulk states appear which displays a uniform density of states throughout the
whole system, see Fig. 2 (d).

Figure 3: $T$-matrix of an impurity at the middle of the lattice with driving
amplitude $A=5.00$. The bulk gap in the spectrum induced by finite spin-orbit
coupling reflects in the $T$-matrix. When the driving frequency is small the
gap remains intact and no states appear within the gap. With increasing
driving frequency states appear in the gap.

The properties of the spectrum can be better understood by considering the
imaginary part of the $T$-matrix, defined by

$\displaystyle{\cal T}=$
$\displaystyle\frac{V_{\text{eff}}}{1-V_{\text{eff}}G^{0}_{0}},$ (13)

as a function of energy. In Fig. 3 we plot $-{\rm Im}\,{\cal T}$ as a function
of the energy for six different frequencies $\Omega$ of the driving field. The
$T$-matrix contains scattering effects of all orders, and the plots in Fig. 3
show $-{\rm Im}\,{\cal T}$ for an impurity located at the center of the
lattice. Since the bulk of a topological insulator does not have any states
available near the Fermi energy, the imaginary part of the $T$-matrix remains
gapped for small driving frequencies. Hence, the spectral content of the
$T$-matrix only marginally deviates from the expected spectral properties of
the pristine lattice. However, huge resonances appear symmetrically outside
the bulk gap, providing coherent resonance peaks. These resonances can be
thought of as in similar terms as the impurity resonances induced in Dirac
materials by particle scattering off local defects AdvPhys.63.1 ;
Nature.403.746 ; PhysRevLett.104.096804 ; PhysRevB.81.233405 ;
PhysRevB.85.121103 ; edge6 ; PhysRevB.94.075401 .

By contrasts, it can be noticed from the plots in Fig. 3 that, with increasing
frequency of the driving field, a finite number of states emerge within the
bulk gap which tends to become filled up by these states. Simultaneously, the
coherent peaks vanish with increased driving frequency as a result of the
necessary charge redistribution which follows from the emergence of additional
resonances in the spectrum.

Next, we consider the impurity to be located at the edge of the lattice, that
is, $(x,y)=(0,25)$, resulting in the plots $-{\rm Im}\,{\cal T}$ shown in Fig.
4 for six different frequencies of the driving field. In stark contrast to the
situation discussed, the coherence peaks outside the bulk gap become are
turned into dips. The metallic nature of the edge states gives rise to finite
scattering states even for small driving frequencies. As the driving frequency
is increased, an increased number of Floquet modes contribute to the
scattering process, which leads to a non-monotonic oscillation pattern in the
$T$-matrix.

Figure 4: $T$-matrix of an impurity at the edge of the lattice with driving
amplitude $A=0.50$. Topological Insulators have edge state at the boundary of
the sample and a finite no of scattering states appear within the bulk gap.
$T$-matrix also, show similar behavior with a finite weight within the gap
even for small frequency region.

##  V Local density of states

In this section, we consider the effect of impurity on the local density of
electron states, something which might be useful for local probing
experiments, such as scanning tunneling microscopy. First, we focus on the
case where the impurity is located at the center of the lattice. In Fig. 5, we
plot the local density of electron states at the impurity site for a sequence
of different values of the amplitude $A$ of the driving field for a fixed
frequency $\Omega$. In the absence of the impurity, $A=0$, the local density
of states at the center (in bulk, far from the edge) of the lattice has a gap
(not shown), which is caused by the finite spin-orbit coupling. For a weak
driving amplitude, the local density of electron states shows a hard gap,
which is reminiscent of the bulk gap in the unperturbed case. Surprisingly,
however, even by increasing the driving amplitude by more than an order of
magnitude, the local density of electron states does not change appreciably
within the gap. This difference is in stark contrast when compared to the
T-matrix results, where a finite weight appears when the impurity is located
at the center of the lattice. The local density of electron states increases
away from the Fermi energy, which can be seen by zooming in inside the gap,
while the local density of electron states remains constant at the Fermi
energy. This implies that even though there is a finite number of states
appearing within the bulk gap due to impurity scattering, the density of
states at the Fermi level does not change. The coherence peaks emerging at the
edges of the gap, become less prominent with increasing amplitude of the
driving field.

Figure 5: Local density of electron states for a single impurity with driving
frequency $\Omega=0.50$ at the impurity site which is located in the middle of
the lattice for small energy. The inset figure shows the local density of
electron states in a large energy range.

In Fig. 6, we plot the local density of electron states for the case where the
impurity is located at the edge of the lattice. It is clear that this set-up
leads to completely different kinds of features in the local density of
electron states, as compared to the previous one. The non-trivial nature of
the metallic TI metallic edge states, due to the chirality, becomes more
significant at the boundary of the lattice. As can be seen in Fig. 6, the
local density of electron states is not gapped in the absence of impurities,
$A=0$. This property is preserved for weak driving fields, $0<A\lesssim 2$,
although it is also clear that both the density of states at the Fermi level
as well as the sharp resonances ($|\omega|\approx 1$) become slightly weaker,
there is, nonetheless, a finite local density of electron states inside the
bulk gap. With even larger amplitudes, the overall low energy local density of
electron states becomes simultaneously broadened and suppressed in amplitude.

##  VI Summary and conclusions

In this paper, we have studied the effects of a time-dependent potential
impurity in two-dimensional topological insulators. By employing a numerically
exact matrix continued fraction method, which takes into account higher order
Floquet modes, we capture the low energy physics throughout a wide range
between weak and strong driving fields. This approach is implemented for both
continuum and lattice models of the topological insulator.

Regarding the continuum model, we find that the low lying energy states remain
unaffected by driving fields with weak to moderate strength, whereas the
higher lying energy states become strongly modified. The linear property of
the local density of electron states disappears with increasing amplitude of
the driving field. Many more oscillations emerge in the local density of
electron states due to contributions from the higher Floquet modes at larger
driving amplitude, as well as a slight increase in the bandwidth.

Figure 6: Local density of electrons states for a single impurity with driving
frequency $\Omega=2.00$ at the impurity site which is located at the edge of
the lattice for small energy. Inset figure show local density of electron
states in a large energy range.

We have imposed open boundary condition for the two-dimensional lattice model
of the topological insulator along with $x$-direction, which gives metallic
edge states. Here, we have analyzed the effect of an impurity at two different
locations, namely, at the center and an edge of the lattice. The bare local
density of electron states has insulating and metallic properties at these two
lattice sites, respectively. The imaginary part of the $T$-matrix has very
different responses for impurities located at these two positions. For a
weakly perturbing impurity located at the center of the lattice, the bulk gap
of the lattice is preserved in the $T$-matrix, while a finite density grows up
within the gap at stronger driving fields caused by increased impurity
scattering. For strong enough fields, this density eventually fills up the
gap. On the other hand for an impurity located at the edge of the lattice, the
$T$-matrix and, hence, the impurity scattering, largely impacts the electronic
density around the Fermi energy already for weak driving fields.

We conclude that despite the apparent simplicity of the model and some
similarities with the physics originating from single impurity interacting
with the surfaces states of a three dimensional topological insulator, there
is increased complexity in the present set-up which can be directly linked to
the time-dependent potential. This should have bearing on issues related to
vibrational defects as well as the external driving field we have in mind
here. However, since the driving frequencies in the two situations should be
expected to differ by orders of magnitude, we can, nonetheless, conjecture
that our low frequency results should be applicable to vibrational impurities.
We anticipate that studies of the time-dependent using different methods
should corroborate our findings.

###### Acknowledgements.

We thank Stiftelsen Olle Engqvist Byggmästare and Vetenskapsrådet for
financial support.

## References

  * (1) D.F. Escande, Phys. Rep 121, 165 (1985)
  * (2) L.E. Reichl and W.M. Zheng, Directions in chaos, 1, Hao Bailin (ed.). Singapore: World Scientific 1987
  * (3) A.J. Lichtenberg and M.A. Lieberman, Regular and stochastic motion. New York: Springer 1981
  * (4) G.M. Zaslavsky, Phys. Rep. 80,157 (1981)
  * (5) A. G. Fainshtein, N. L. Manakov, V. D. Ovsiannikov, and L. P. Rapoport, Phys. Rep. 210, 111 (1992).
  * (6) F. L. Moore, J. C. Robinson, C. Bharucha, P. E. Williams, and M. G. Raizen, Phys. Rev. Lett. 73, 2974 (1994);
  * (7) P. J. Bardroff, I. Bialynicki-Birula, D. S. Krähmer, G. Kurizki, E. Mayr, P. Stifter, and W. P. Schleich, Phys. Rev. Lett. 74 3959 (1995).
  * (8) N. H. Lindner, G. Refael, and V. Galitski, Nat. Phys. 7, 490(2011).
  * (9) T. Oka and H. Aoki, Phys. Rev. B 79, 081406 (2009).
  * (10) T. Kitagawa, T. Oka, A. Brataas, L. Fu, and E. Demler, Phys. Rev. B 84, 235108 (2011).
  * (11) L. Jiang, T. Kitagawa, J. Alicea, A. R. Akhmerov, D. Pekker, G. Refael, J. I. Cirac, E. Demler, M. D. Lukin, and P. Zoller, Phys. Rev. Lett. 106, 220402 (2011).
  * (12) N. H. Lindner, D. L. Bergman, G. Refael, and V. Galitski, Phys. Rev. B 87, 235131 (2013).
  * (13) Y. T. Katan and D. Podolsky, Phys. Rev. Lett. 110, 016802 (2013).
  * (14) G. Usaj, P. M. Perez-Piskunow, L. E. F. Foa Torres, and C. A. Balseiro, Phys. Rev. B 90, 115423 (2014).
  * (15) Z. Gu, H. A. Fertig, D. P. Arovas, and A. Auerbach, Phys. Rev. Lett. 107, 216601 (2011).
  * (16) D. E. Liu, A. Levchenko, and H. U. Baranger, Phys. Rev. Lett. 111, 047002 (2013).
  * (17) Y. Li, A. Kundu, F. Zhong, and B. Seradjeh, Phys. Rev. B 90, 121401 (2014).
  * (18) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005). (2005).
  * (19) J. H. Shirley, Phys. Rev. 138, B979 A965).
  * (20) Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik, Science 342, 453 (2013).
  * (21) F. Mahmood, C.-K. Chan, Z. Alpichshev, D. Gardner, Y. Lee, P. A. Lee, and N. Gedik, Nat. Phys. 12, 306 (2016).
  * (22) M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, Nature (London) 496, 196 (2013).
  * (23) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Nature (London) 515, 237 (2014).
  * (24) H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nat. Phys. 5, 438 (2009).
  * (25) Y. L. Chen et al., Science 325, 178 (2009).
  * (26) Y. Xia et al., Nat. Phys. 5, 398 (2009). (2007).
  * (27) W.-C. Lee, C. Wu, D. P. Arovas, and S.-C. Zhang, Phys. Rev. B 80,245439 (2009);
  * (28) X. Zhou, C. Fang, W.-F. Tsai, and J. P. Hu, Phys. Rev. B 80, 245317 (2009);
  * (29) H.-M. Guo and M. Franz, Phys. Rev. B 81, 041102 (2010);
  * (30) A. M. Black-Schaffer and A. V. Balatsky, Phys. Rev. B 85, 121103(R) (2012);
  * (31) J. Fransson, A. M. Black-Schaffer, and A. V. Balatsky, Phys. Rev. B, 90, 241409(R) (2014);
  * (32) A. M. Black-Schaffer, A. V. Balatsky, and J. Fransson, Phys. Rev. B 91, 201411(R) (2015).
  * (33) J. Fransson, A. M. Black-Schaffer, and A. V. Balatsky Phys. Rev. B 94, 075401 (2016).
  * (34) P. Roushan, J. Seo, C. V. Parker, Y. S. Hor, D. Hsieh, D. Qian, A. Richardella, M. Z. Hasan, R. J. Cava, and A. Yazdani, Nature (London) 460 (2009);
  * (35) T. Zhang, P. Cheng, X. Chen, J.-F. Jia, X. Ma, K. He, L. Wang, H. Zhang, X. Dai, Z. Fang et al., Phys. Rev. Lett. 103, 266803 (2009);
  * (36) Z. Alpichshev, J. G. Analytis, J.-H. Chu, I. R. Fisher, Y. L. Chen, Z. X. Shen, A. Fang, and A. Kapitulnik, Phys. Rev. Lett. 104, 016401 (2010).
  * (37) A. Roth, C. Brüne, H. Buhmann, L. W. Molenkamp, J. Maciejko, X.-L. Qi, and S.-C. Zhang, Science 325 , 294 (2009);
  * (38) M. König, S. Wiedmann, C. Brüe, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007);
  * (39) G. M. Gusev, Z. D. Kvon, O. A. Shegai, N. N. Mikhailov, S. A. Dvoretsky, and J. C. Portal, Phys. Rev. B 84, 121302(R) (2011).
  * (40) G. Casati, I. Guarneri, D. L. Shepelyansky, IEEE J. Quant. Elect. 24, 1420 (1988);
  * (41) R. Shakeshaft, Comm. At. Mol. Phys. 28, 179 (1992);
  * (42) Laser Phys. 3, No. 2 (1993), special issue on Atoms, Ions and Molecules in a Strong Laser Field, edited by A. M. Prokhorov.
  * (43) G. Floquet, Ann. de l’Ecole Norm. Sup. 12, 47 (1883);
  * (44) E. L. Ince, Ordinary Differential Equations (Dover Publ., New York, 1956);
  * (45) W. Magnus and S. Winkler, Hill’s Equation (Dover Publ., New York, 1979).
  * (46) C. Wu, B. A. Bernevig, and S.-C. Zhang, Phys. Rev. Lett. 96,106401(2006).
  * (47) C. Xu and J. E. Moore, Phys. Rev. B 73, 045322 (2006).
  * (48) Risken. H, The Fokker-Planck equation. Springer Series in Synergetics, Vol. 18. Berlin, Heidelberg, New York: Springer 1984
  * (49) D F Martinez, J. Phys. A: Math. Gen. 36 9827 (2003).
  * (50) F. Grossmann, P. Jung, T. Dittrich, and P. Hänggi, Z. Phys. B 84, 315 (1991).
  * (51) J. -H. She, J. Fransson, A. R. Bishop, and A. V. Balatsky, Phys. Rev. Lett. 110, 026802 (2013).
  * (52) J. Fransson, J -H. She, L. Pietronero, and A. V. Balatsky, Phys. Rev. B, 87, 245404 (2013).
  * (53) T. Wehling, A. Black-Schaffer, and A. Balatsky, Adv. Phys. 63, 1 (2014).
  * (54) S. H. Pan, E. W. Hudson, K. M. Lang, H. Eisaki, S. Uchida, and J. C. Davis, Nature (London) 403, 746 (2000).
  * (55) M. M. Ugeda, I. Brihuega, F. Guinea, and J. M. Gómez-Rodríguez, Phys. Rev. Lett. 104, 096804 (2010).
  * (56) R. R. Biswas and A. V. Balatsky, Phys. Rev. B 81, 233405 (2010).
  * (57) A. M. Black-Schaffer and A. V. Balatsky, Phys. Rev. B 85, 121103 (2012).


11institutetext: 1: Department of Astronomy, University of Michigan, Ann
Arbor, MI, USA

# On The Orbital Separation Distribution
and Binary Fraction of M Dwarfs

Nicholas Susemiehl 11 Michael R. Meyer 11

(Accepted September 9, 2021)

###### Abstract

Aims. We present a new estimate for the binary fraction (the fraction of stars
with a single companion) for M dwarfs using a log-normal fit to the orbital
separation distribution.

Methods. We use point estimates of the binary fraction (binary fractions over
specific separation and companion mass ratio ranges) from four M dwarf surveys
sampling distinct orbital radii to fit a log-normal function to the orbital
separation distribution. This model, alongside the companion mass ratio
distribution given by Reggiani & Meyer (2013), is used to calculate the
frequency of companions over the ranges of mass ratio (q) and orbital
separation (a) over which the referenced surveys were collectively sensitive -
[0.60 $\leq$ q $\leq$ 1.00] and [0.00 $\leq$ a $\leq$ 10,000 AU]. This method
was then extrapolated to calculate a binary fraction which encompasses the
broader ranges of [0.10 $\leq$ q $\leq$ 1.00] and [0.00 $\leq$ a ¡ $\infty$
AU]. Finally, the results of these calculations were compared to the binary
fractions of other spectral types.

Results. The binary fraction over the constrained regions of [0.60 $\leq$ q
$\leq$ 1.00] and [0.00 $\leq$ a $\leq$ 10,000 AU] was calculated to be
$0.229\pm 0.028$. This quantity was then extrapolated over the broader ranges
of q (0.10 - 1.00) and a (0.00 - $\infty$ AU) and found to be
$0.462^{+0.057}_{-0.052}$. We used a conversion factor to estimate the
multiplicity fraction from the binary fraction and found the multiplicity
fraction over the narrow region of [0.60 $\leq$ q $\leq$ 1.00] and [0.00
$\leq$ a $\leq$ 10,000 AU] to be $0.270\pm 0.111$. Lastly, we estimate the
multiplicity fractions of FGK, and A stars using the same method (taken over
[0.60 $\leq$ q $\leq$ 1.00] and [0.00 $\leq$ a $\leq$ 10,000 AU]) and find
that the multiplicity fractions of M, FGK, and A stars, when considered over
common ranges of q and a, are more similar than generally assumed.

###### Key Words.:

(stars:) binaries: general – stars: late-type – stars: statistics

##  1 Introduction

Because M dwarfs comprise as much as 75% of local stellar populations (Henry
et al., 2006; Winters et al., 2019, and references therein), it is important
to have an accurate census of their properties. Understanding the multiplicity
fraction, the fraction of stellar systems with one or more gravitationally-
bound companion (regardless of whether it is above or below the hydrogen
burning limit), of M dwarfs will have significant implications for star
formation theories, inform us about the emergence of planetary systems around
low mass stars, and allow for the modeling of both galactic and extra-galactic
stellar populations. A related quantity is the binary fraction, the fraction
of stars with a single gravitationally-bound companions. The key difference
between these two statistics is that the multiplicity fraction considers
systems with more than two companions while the binary fraction does not.
Mathematically, these quantities are expressed as
$MF=\frac{N_{D}+N_{T}+N_{H}}{N_{S}+N_{D}+N_{T}+N_{H}}$ and
$BF=\frac{N_{D}}{N_{S}+N_{D}+N_{T}+N_{H}}$ (where $MF$ denotes the
multiplicity fraction, $BF$ the binary fraction, and $N_{S}$, $N_{D}$,
$N_{T}$, $N_{H}$ the number of single, double, triple, and higher-order
systems respectively). These differences extend to formation theories: the
multiplicity fraction may be more related to gravitational fragmentation of
star forming molecular cloud cores than the more simple binary fraction (cf.
Raghavan et al., 2010; De Rosa et al., 2014). In addition, differences between
the multiplicity and binary fractions as a function of primary star mass are
another key constraint of formation theory. Multiple conventions exist to
define the mass ratio and orbital distribution for higher order systems so we
only consider the binary fraction when fitting our model in order to avoid
these complications. Despite its import, the multiplicity properties of M
dwarfs are still not yet fully understood.

We model the binary fraction using the companion mass ratio distribution
($\psi$) and orbital separation distribution ($\phi$). The companion mass
ratio distribution represents the fraction of stars which have a companion
within a given range of the mass ratio. The mass ratio, $q$, is defined as the
ratio between the mass of a secondary star and the mass of the respective
primary star where, by definition, the mass of the primary star is greater
than that of the secondary so that $q\leq 1$. The orbital separation
distribution is the fraction of stars which have a companion as a function of
the separation between the two stars. The semi-major axis of an orbital
ellipse, $a$, serves as a measure of the separation between the primary and
secondary stars. Integrating either of these functions over specific ranges of
their respective parameters returns the fraction of M dwarfs which have a
single companion over that range if one assumes the other parameter has been
integrated over its full parameter space. $\psi$ and $\phi$ must be combined
to provide a binary fraction over a fixed ranges of q and a, calculated as:

$f=\int_{q_{min}}^{q_{max}}\psi
dq*\int_{log_{10}a_{min}}^{log_{10}a_{max}}\phi dlog_{10}(a)$ (1)

where $q_{min},q_{max},log_{10}a_{min},$ and $log_{10}a_{max}$ represent the
lower and upper bounds for the mass ratio and semi-major axes of interest.

The companion mass ratio distribution is discussed in Reggiani & Meyer (2013),
who find:

$\psi=q^{.25\pm 0.29}$ (2)

describes this distribution for M dwarfs as well as field FGK stars and
intermediate mass stars in the young ScoOB2 association. The present work
attempts to fit the functional form $\phi$ for the orbital separation
distribution.

This method for calculating the binary fraction is built upon a key
assumption: that the companion mass ratio distribution does not depend on
orbital separation. This assumption allows for the calculation of binary
fractions using Equation (1). Reggiani & Meyer (2013) investigate this
assumption for M dwarfs based on the data of Janson et al. (2012) and are
unable to reject it. Furthermore, we only take into consideration binary
systems and not triplets, quadruplets, etc. We further discuss these caveats
in Section 4.

Numerous attempts have been made to constrain the multiplicity of M dwarfs
(e.g. Fischer & Marcy, 1992; Janson et al., 2012; Winters et al., 2019), but
none include samples of stellar binary pairs that are complete over the full
ranges of mass ratio ($0\leq q\leq 1$) or separation ($0<a<\infty$). Instead,
observational xresults are necessarily limited to specific ranges of
parameters. These past works have made useful constraints to important
parameters regarding the population of M dwarfs. For instance, Fischer & Marcy
(1992) found that about $42\pm 9\%$ of nearby M dwarfs (with mass ratio
0.2¡q¡1 and separation 0 ¡a¡ 10,000 AU) have a companion and Ward-Duong et al.
(2015) estimated a companion star fraction of $23.5\%\pm 3.2\%$ over a
similarly wide region of the parameter space (3¡a¡10,000 AU, 0.2¡q¡1). Janson
et al. (2012) suggests that the multiplicity of M dwarfs is intermediate
between that of higher-mass stars and brown dwarfs. Cortés-Contreras et al.
(2017) found a peak in the orbital separation distribution of M dwarf multiple
system at 2.5 - 7.5 AU, which is similar to the broad peak at around 4-20 AU
found by Winters et al. (2019) (cf. Winters et al., 2021).

In Section 2 of this paper, we describe the methods used to fit the orbital
separation distribution and calculate the binary fraction. In Section 3, the
results are presented. In Section 4, we discuss the implications and further
investigate the assumptions used. Section 5 summarizes our work.

##  2 Methods

###  2.1 The data

The first step in exploring the orbital separation distribution and binary
fraction of M dwarfs was to synthesize point estimates of the binary fraction
from a variety of different M dwarf surveys. Data was compiled from surveys
which employed the radial velocity (Delfosse et al., 1999) and direct imaging
(Cortés-Contreras et al., 2017; Janson et al., 2012; Ward-Duong et al., 2015)
companion detection methods. This allowed us to investigate companions over a
broad range of semi-major axes. The projected separations from the direct
imaging surveys of Cortés-Contreras et al. (2017) and Ward-Duong et al. (2015)
were converted to physical separations using the correction factor described
in Fischer & Marcy (1992) of 1.26 (Janson et. al. 2012 already made a
correction using this factor). Other works, including Dupuy & Liu (2011),
suggest alternative forms to this correction factor that account for other
orbital parameters such as eccentricity, but the correction factor estimated
by Dupuy & Liu (2011) ($1.16^{+0.81}_{-0.31}$) is similar to the 1.26
suggested by Fischer & Marcy (1992).

We include only detected multiple systems associated with values of mass ratio
and semi-major axis over which each survey was at least 90% complete. The
lower limit on mass ratio considered was set by Cortés-Contreras et al. (2017)
as q $\geq$ 0.60. Any detected multiple system which was separated by a
distance outside the range that each particular survey was at least 90%
complete was also removed (i.e. no longer considered a detection of a binary
system). We also excluded any higher order systems by dropping such systems
from what we consider to be companions but including them in the total primary
star number, effectively treating these systems as single stars. Delfosse et
al. (1999) detected six higher order systems, Cortés-Contreras et al. (2017)
detected four, Janson et al. (2012) detected 14, and Ward-Duong et al. (2015)
detected five. The binary fraction point estimates for each survey were
calculated by dividing the remaining number of binary pairs by the total
parent population of each survey. The error bars on the binary fraction point
estimates were taken to be the Poisson counting error of the respective
surveys (calculated after the aforementioned cleaning steps). Table (1) gives
the point estimate of the binary fraction for each survey along with the
respective detection method and range of semi-major axis. Similarly, Table (2)
depicts the total number of parent stars surveyed and the number of companion
detections for each reference.

In certain cases, assumptions had to be made about details of the surveys.
Delfosse et al. (1999) reported the period and orbital velocity of stars in
their sample, but not the semi-major axis. Kepler’s Third Law and the masses
of the primaries were used to convert the given values into semi-major axes.
Janson et al. (2012) did not note the minimum mass ratio that the survey could
detect. The authors state: ”In every case where $q_{m}*m_{a}<0.08$ $M_{sun}$,
the system is removed from the analysis… Furthermore, we include only stars
with primary mass ¿ 0.2 $M_{sun}$, since lower mass stars are not fully
complete out to 52 pc” (Janson et al., 2012). From this, we took a minimum
secondary mass of 0.08 $M_{sun}$ and a minimum primary mass of 0.2 $M_{sun}$
to calculate a minimum mass ratio of 0.4. Both Janson et al. (2012) and Ward-
Duong et al. (2015) included a small number of multiple systems with K-type
primary stars, each of which was removed from this analysis. The data from
Ward-Duong et al. (2015) was split into two bins by orbital separation: 0 -
100 AU and 100 - 10,000 AU because it covers such a wide range of semi-major
axes (3 - 10,000 AU).

In order to justify our use of the Reggiani & Meyer (2013) power law (Equation
2), we employed the Kolmogorov–Smirnov (KS) test to compare the mass ratio
samples of Ward-Duong et al. (2015) and Cortes-Contreras et al. (2016) to ones
drawn from the mass ratio power law distribution of Reggiani & Meyer (2013)
(considering only q¿0.6; Fischer & Marcy 1992 and Janson et al. 2012 were used
to fit the power law of Reggiani and Meyer 2013 and so were not tested here).
These KS tests returned p-values of 0.04 and 0.02 respectively. Based on these
p-values, we do not reject the null hypothesis that there is no difference
between the distributions of Ward-Duong et al. (2015) and Cortes-Contreras et
al. (2016) and the mass ratio distribution proposed by Reggiani & Meyer (2013)
at the 0.01 significance level. As an alternative, future similar work could
consider using the mass ratio distribution offered by El-Badry et al. (2019)
which suggests a broken power law fit and found an excess of near-equal mass
systems. Using this mass ratio distribution would likely change the results
found in this study.

Table 1: Binary Fraction Point Estimates, q $\geq$ 0.6

Reference | Semi-Major Axis Range (AU) | Primary Mass Range ($M_{sun}$) | Binary Fraction Estimate
---|---|---|---
Delfosse et al. (1999) | 0.00 - 4.63 | 0.13 - 0.44 | 0.04 $\pm$ 0.018
Cortés-Contreras et al. (2017) | 3.28 - 37.17* | 0.14 - 0.49 | 0.07 $\pm$ 0.012
Janson et al. (2012) | 3.00 - 227.00 | 0.20 - 0.59 | 0.17 $\pm$ 0.015
Ward-Duong et al. (2015) A | 3.78 - 126.00* | 0.24 - 0.67 | 0.11 $\pm$ 0.022
Ward-Duong et al. (2015) B | 126.00 - 12,600.00* | 0.32 - 0.52 | 0.07 $\pm$ 0.017

111This table depicts each of the M dwarf surveys used in this study alongside
their respective ranges of semi-major axis over which the survey was at least
90% complete, ranges of primary star masses, and the point estimates of the
binary fraction that we made after excluding companion detections outside the
respective semi-major axis ranges and with q ¡ 0.6.
* Projected separations converted into physical separations by multiplying limits by 1.26 (Fischer & Marcy, 1992).

Table 2: Survey Details

Reference | Total Num. of Stars Surveyed | Total Num. of Detected Companions | Companion Subsample
---|---|---|---
Delfosse et al. (1999) | 127 | 12 | 5
Cortés-Contreras et al. (2017)* | 490 | 80 | 34
Janson et al. (2012) | 701 | 219 | 117
Ward-Duong et al. (2015) A | 232 | 65 | 25
Ward-Duong et al. (2015) B | 232 | 65 | 16

222This table summarizes each of the M dwarf surveys used in this study. The
first column lists the total number of type M primary stars included in each
survey. The second column is the number of companions detected by each survey.
The third column is the number of companions remaining after the corrections
described in Section 2.1 are made. The binary fraction estimates described in
Table 1 can be found by dividing the fourth column of this table by the
second.
*:We used the volume-limited sub-sample of Cortés-Contreras et. al. (2017).

###  2.2 Fitting the model

We fit a log-normal model to the orbital separation distribution of the M
dwarf binary systems by using the binary fraction point estimates from the
referenced surveys (Table 1) as benchmarks to those computed using Equation
(1). This model (Equation 3) has as free parameters an amplitude (A) which
normalizes the function, a mean ($\mu$) which indicates the peak of the
distribution, and a standard deviation ($\sigma$) which describes the spread
of the distribution.

$\phi=A*\frac{e^{-(log_{10}(a)-\mu)^{2}/(2\sigma^{2})}}{\sigma*\sqrt{2\pi}}$
(3)

In order to fit the model and find the best values for these parameters, a
Markov Chain Monte Carlo (MCMC) routine was utilized via the emcee Python
package (Foreman-Mackey et al., 2013). Model binary fractions, which are
calculated via Equation (1) by inputting parameter values into $\phi$, were
compared to the binary fraction point estimates from the referenced surveys
using a chi-squared log-likelihood function. Uniform log-priors were chosen
with bound $0.0\leq A\leq 2.0$, $-1.0\leq\mu\leq 2.0$, and $-1.0\leq\sigma\leq
2.0$ (the latter two given in the log of the semi-major axis in AU). The MCMC
was initialized with 100 walkers positioned randomly within the bound of the
log-prior and run for 5,000 steps. After this, the flattened chains of each of
the three parameters were plotted in order to visualize the convergence of the
walkers. This informed the decision to discard the first 200 steps of each
flatchain as the burn-in values. The best fit amplitude, mean, and standard
deviation were then calculated as the median value of these chains, with error
corresponding to the 68% confidence interval of the respective posterior
distributions. A cornerplot (Figure 1) was produced which shows the
distributions of each of the three parameters.

Figure 1: Marginal distributions of each of the three parameters as well as
the covariances between all of them. The median values and 68% confidence
intervals of each of the parameters are drawn as vertical lines and also noted
above the respective histograms.

###  2.3 Calculating the binary fraction

Next, the model for the orbital separation distribution, with its best-fitting
parameters, was used alongside the companion mass ratio distribution model
from Reggiani & Meyer (2013) to calculate the binary fraction of M dwarfs as:

$\begin{split}f=\int_{q_{min}}^{q_{max}}q^{.25}dq*\\\
A_{best}*\int_{log_{10}a_{min}}^{log_{10}a_{max}}\frac{e^{-(log_{10}(a)-\mu_{best})^{2}/(2\sigma_{best}^{2})}}{\sigma_{best}*\sqrt{2\pi}}dlog_{10}(a)\end{split}$
(4)

This calculation resulted in a binary fraction that is representative over
these limited ranges of mass ratio and semi-major axis. These ranges were
later expanded to [0.10 $\leq$ q $\leq$ 1.0] and [0 $\leq$ a ¡ $\infty$ AU] to
allow for the calculation of a broader M dwarf binary fraction. The error on
the binary fraction was calculated as the 68% confidence interval of the
probability distribution function of the binary fraction, which was generated
by sampling 1,000 values of A, $\mu$, and $\sigma$ from the posterior space
and calculating a value of the binary fraction for each of these parameter
values (as per Equation 4). Figure (2) summarizes the model’s fit to the data.

Figure 2: Comparison of the binary fraction survey data to the fitted model.
Each of the red horizontal lines represent one of the five binary fraction
estimates from the literature. The right vertical axis position of these lines
is the binary fraction reported in the data and the span along the horizontal
axis of these lines is the range in semi-major axes that the respective data
cover. The red vertical lines are the error in these binary fraction
estimates. The blue curve is the fitted log-normal model of the orbital
separation distribution. Each red dot is the binary fraction estimate obtained
by integrating the fitted model over the range of semi-major axis respective
to each referenced survey. When visually assessing the validity of our model,
the most important aspect of this plot to notice is that each red data point
fits within each red vertical line, implying that our model replicates the
referenced binary fraction point estimates within error.

##  3 Results

###  3.1 Orbital separation model

The best-fit to point estimates of the binary fraction from the four M dwarf
surveys resulted in the parameters of A, $log_{10}(\mu)$, and
$log_{10}(\sigma)$, being $0.61\pm 0.07$, $1.68^{+0.14}_{-0.16}$, and $0.97\pm
0.19$ respectively. This fit is associated with a reduced chi-squared
parameter value of 1.315. With two degrees of freedom, the probability of
achieving a value greater than or equal to 1.315 is 0.269. Therefore, we could
not reject the null hypothesis that the data came from this model.

###  3.2 Binary fraction

Integrating Equation (4) over the constrained regions of [0.60 $\leq$ q $\leq$
1.00] and [0 $\leq$ a $\leq$ 10,000 AU] resulted in a binary fraction of
$0.229\pm 0.028$. Equation (4) was then integrated again over [0.1 $\leq$ q
$\leq$ 1.0] and [0.0 $\leq$ a ¡ $\infty$ AU] to calculate a broad binary
fraction of $0.462^{+0.057}_{-0.052}$. This binary fraction considers stellar
and sub-stellar companions to M dwarfs, but not planetary companions. The
binary fraction corresponding to mass ratios and semi-major axes outside of
the ranges of [0.60 $\leq$ q $\leq$ 1.00] and [0 $\leq$ a $\leq$ 10,000 AU]
must be considered with care because they are extrapolations of our model.

###  3.3 Multiplicity fraction

Due to its import in constraining theories of stellar formation (Janson et
al., 2012; De Rosa et al., 2014; Winters et al., 2019, and others), we provide
an estimate of the multiplicity fraction based on our binary fraction fit.
Winters et al. (2019) estimates the ratio of multiple systems (including
binary and higher order systems) to binaries to be 1.179 $\pm$ 0.107 (with
error propagated from the Poisson error in the number of multiple and binary
systems). By assuming that the binary fraction and the total number of higher
order systems are independent, we were able to multiply our binary fraction by
this ratio to obtain a multiplicity fraction estimate of $0.270\pm 0.111$ over
the constrained regions of [0.60 $\leq$ q $\leq$ 1.00] and [0 $\leq$ a $\leq$
10,000 AU]. The multiplicity of M dwarfs integrated over the same region of
mass ratio but over all values of separation (0-$\infty$) is presented here
for use in the comparison made in section 4.2 to be 0.272 $\pm$ 0.111.

##  4 Discussion

###  4.1 Comparisons to other M dwarf surveys

We compare our model to the to the observation of Bowler et al. (2015), which
estimated the frequency of brown dwarf companions to M dwarf hosts. They
observe the frequency of brown dwarf companions to M dwarf primaries over a
semi-major axis range of 10 - 100 AU and mass ratio range of 0.039 - 0.224 to
be $0.028^{+0.024}_{-0.015}$. We integrated Equation (1) with the best fit
parameters described above over these ranges of semi-major axis and mass ratio
and obtained an binary fraction estimate of $0.027^{+0.010}_{-0.007}$ (cf.
Dieterich et al., 2012). While the sample of Bowler et al. (2015) is drawn
from nearby young moving groups and our preceding work was built using samples
drawn from field stars, this comparison is justified because the multiplicity
of stars in these moving groups does not appear to be discernibly different
from the field (see Elliott et al., 2015; De Furio et al., 2019). None of the
surveys used to construct our binary fraction model probed sub-stellar
companions, and so the capability of our model to reproduce the results of
Bowler et al. (2015) is of particular note.

Winters et al. (2019) (W19) also studied the orbital separation distribution
and binary fraction of M dwarfs. They find an upper limit to the peak in the
this distribution to be at 20 AU (projected separation) and note that most
stellar companions to M dwarfs lay between 4 and 20 AU. When corrected to a
physical separation using the 1.26 correction factor of Fischer & Marcy (1992)
(Section 2.1), W19’s peak to the orbital separation distribution becomes 25
AU. This falls about 1.9 standard deviations from our peak estimate of
$48^{+18}_{-15}$ AU (our orbital separation model estimates the probability of
the peak of this distribution occurring at less than 20 AU to be about 3%.) An
attempt was made to draw a comparison between our binary fraction and that of
W19. However, it was not straight-forward to find the range of q and
separation over which the survey was complete. We therefore attempted to
estimate this ourselves. First, we calculated a binary fraction from the W19
sample by using their stated number of single:double:triple:high-order
systems. We assumed that the binary fraction uniformly covered mass ratios
from a minimum q set by the hydrogen burning limit divided by the median
primary mass of their sample to a maximum q of 1. We assumed the companion
mass ratio distribution of Reggiani & Meyer (2013) to convert this to the
range of mass ratio that is representative of our binary fraction, 0.6 - 1.0.
Furthermore, we assumed that the multiplicity fraction of W19 was
representative for companions with separations less than 10,000 AU. However,
the binary fraction that we estimated from the sample of W19, 0.133 $\pm$
0.014 did not fall within error of our constrained binary fraction, 0.229
$\pm$ 0.028. This discrepancy is likely due to the numerous assumptions we
were forced to make.

Another survey of M dwarf multiple systems, Shan et al. (2015), finds the
fraction of M dwarf companions with periods less than 90 days to M dwarf hosts
to be $0.11^{+0.02}_{-0.04}$. This estimate is much higher than Delfosse et
al. for closely-separated M dwarfs and would have a significant impact on our
fit if it were included as a data point. However, we were unable to include
Shan et al. (2015) in our analysis because the orbital separation and mass
ratio completeness limits and values for each system are not easy to
ascertain. Jódar et al. (2013) was also considered for inclusion in our fit
but it was also not clear over what range of mass ratios and orbital
separations to which the survey was complete. The microlensing survey
Shvartzvald et al. (2016) was not used in this analysis because there were
difficulties in translating this work into usable point estimates without
making strong assumptions that were not able to verify regarding their
sensitivity limits in mass ratio and orbital separation. The inclusion of
microlensing surveys would improve our analysis significantly because of their
sensitivity to low-separation companions

###  4.2 Comparisons to other spectral type binary fractions

We now compare the binary fraction of M dwarfs to that of sun-like FGK- and
more massive A-type stars over the range of mass ratios studied here but over
all separations (0.6 $\leq$ q $\leq$ 1.0 and 0.00 $\leq$ a $\leq$ $\infty$
AU). We expand this integration to be over all separations because of the
greater prevalence of higher mass stars at larger separations. In order make
this comparison, we choose an appropriate form of the companion mass ratio
distribution for the stellar hosts in question then recalculate the orbital
separation distribution for each population of stars based on previous works.
To compare to the FGK binary fraction, we extrapolated the model of orbital
separation distribution from Raghavan et al. (2010) (which peaks at 48 AU) to
integrate it over separations from 0 to $\infty$ AU. This was used alongside
the companion mass ratio distribution of Reggiani and Meyer (2013) integrated
from 0.6 to 1.0 to find an FGK-star binary fraction of 0.17 $\pm$ 0.03, where
the error is estimated as the Poisson counting error based on the survey of
Raghavan et al. (2010). Following the process in Section 3.3 of this work, we
convert this binary fraction to a multiplicity fraction using the ratio of
multiple systems to binary systems described in Raghavan et al. (2010) (1.33
$\pm$ 0.14) and find an FGK-star multiplicity fraction of $0.23\pm 0.14$ over
0.6 $\leq$ q $\leq$ 1.0 and 0.00 $\leq$ a $\leq$ $\infty$ AU. We made a
similar calculation for the A star binary fraction, this time based on De Rosa
et al. (2014) the orbital separation distribution of this survey has a mean at
nearly 400 AU). However, because the A star sample of Reggiani & Meyer (2013)
(young stars in an association) is different than that of De Rosa et al.
(2014) (field stars) we chose to instead use the companion mass ratio
distribution distribution described in De Rosa et al. (2014) to make this
calculation (which is similar to that of Reggiani & Meyer (2013) at close
separations but not at separations greater than 125 AU). Using this, we found
an A star binary fraction of $0.11\pm 0.02$ and multiplicity fraction of
$0.15^{+0.25}_{-0.15}$ (by multiplying the binary fraction by the ratio of
multiple systems to binary systems of De Rosa et al. (2014), 1.36 $\pm$ 0.25)
over mass ratio 0.6 $\leq$ q $\leq$ 1.0 and orbital separation 0.00 $\leq$ a
$\leq$ $\infty$ AU. The binary fractions calculated here decrease as spectral
type increases, but the multiplicity fractions remain consistent (with large
error bars). These calculations are built on a number of assumptions. Most
notably, that the orbital separation and companion mass ratio distributions
used are complete over the entire space of 0.6 $\leq$ q $\leq$ 1.0 and 0.00
$\leq$ a $\leq$ $\infty$ AU (the orbital separation distribution of De Rosa et
al. (2014) is incomplete at low separations). Also that the multiplicity
fraction can be estimated from the binary fraction as we have done above.
Importantly, these results must be considered with these assumptions and with
their large error bars in mind. This finding is in contrast to Duchêne & Kraus
(2013), which suggests that the occurrence rate of companions increases with
spectral type. Nonetheless, when compared over common ranges of mass ratios,
the multiplicity fraction as a function of spectral type is more similar than
commonly appreciated.

###  4.3 Additional caveats

The binary fraction study conducted here was only concerned with the presence
of binary systems within the M dwarf population. No considerations were made
to account for the presence of higher order (triple, quadruple, etc.) systems
during the fitting process. Instead, we converted the binary fraction to a
multiplicity fraction using the empirical findings of other works. While such
systems do exist, they are relatively rare. A recent M dwarf multiplicity
study found that only 3.3% of M dwarfs exist as triple and higher order
systems (Winters et al., 2019). However, it appears that higher order systems
are more common around higher mass stars (Raghavan et al., 2010; Ward-Duong et
al., 2015).

Our work assumes that the mass ratio of a stellar binary system does not
depend on the separation between its components (Reggiani & Meyer, 2013). Moe
& Di Stefano (2017) found that the distributions of mass ratio and period
(which is directly proportional to separation) of O- and B-type main sequence
primaries are not independent. If this is also true for M dwarf binaries, then
Equation (1) would not hold. For the time being, we accept the findings of
Reggiani & Meyer (2013) but note that further work exploring the
interdependence of mass ratio and orbital separation for specifically type M
stars is needed.

Originally, our fit to the orbital separation distribution utilized the
multiplicity estimate from the radial velocity data of Fischer & Marcy (1992).
The inclusion of this data point (binary fraction of 0.083 $\pm$ 0.03 over 0.6
$\leq$ q $\leq$ 1 and 0.04 $\leq$ a $\leq$ 4) gave a worse fit and led to a
higher reduced chi-squared of 2.143. We believe this is because the primary
mass range of the sample of Fischer & Marcy (1992), 0.14 - 0.33 $M_{Sun}$,
does not include higher mass M dwarfs present in the other referenced surveys.
While this primary mass range may not seem substantially different than that
of Delfosse et al. (1999), we examined the spectral type distributions of the
parent populations of Fischer & Marcy (1992) and Delfosse et al. (1999),
(Marcy & Benitz, 1989; Delfosse et al., 1998, respectively), and found that
they are distributed significantly differently (i.e. Marcy & Benitz 1989, and
therefore Fischer & Marcy 1992, include more low mass stars than Delfosse et
al. 1998/Delfosse et al. 1999) by performing a KS test and receiving a p-value
of 0.0007. Previous work suggests the companion mass ratio distribution for
brown dwarfs is less flat than that of M dwarfs (Burgasser et al., 2007;
Fontanive et al., 2018). The inclusion of Fischer & Marcy (1992) may have been
problematic because its sample is closer to the brown dwarf regime than the
other referenced surveys. Therefore, we chose to exclude this point from our
model fit.

##  5 Summary

We find that the distribution of M dwarf binary separations peaks between
33-66 AU within a 68 % confidence interval (23-85 AU with 95 % confidence),
which is comparable to the peak for FGK stars but smaller than the peak for A
stars (Raghavan et al., 2010; De Rosa et al., 2014). Using the model built
from this finding as well as the companion mass ratio distribution given by
Reggiani & Meyer (2013), we find that the binary fraction over [0.60 $\leq$ q
$\leq$ 1.00] and [0.00 $\leq$ a $\leq$ 10,000 AU] for M dwarfs to be $0.229\pm
0.028$ (this is the most robust of our findings as it required the fewest
extrapolations). We used this to derive the multiplicity fractions over the
same range of mass ratio and over all separations for M, FGK, and A stars of
$0.27\pm 0.11$, $0.23\pm 0.14$, and $0.15\pm 0.25$ respectively. We note that
these three values suggest that the the multiplicity fraction is similar over
these three spectral types and (within these ranges of mass ratio and orbital
separation) but the error bars on these values are large. Furthermore, we find
the binary fraction of M dwarfs over 0.10 $\leq$ q $\leq$ 1.00 and 0.00 $\leq$
a $\leq$ $\infty$ AU to be $0.462^{+0.057}_{-0.052}$, implying that nearly
half of all M dwarfs (over all q and a) have a stellar or sub-stellar
companion. This assertion is highly dependent on the validity of the
assumptions we made. Namely, that the frequency of binary companions is
consistent beyond the parameter space of the surveys used to build our model.
Data which would help improve the predictive power of our orbital separation
model may include M dwarf multiplicity surveys which probe pairs with
separations between 3 and 30 AU and low mass ratios. Future studies may expand
upon this work by surveying the binary fraction of M dwarf companions in the
brown dwarf regime, explore correlations between the companion mass ratio and
orbital separation distributions, as well as compare M dwarf binary fraction
to L, T, and Y- Dwarf multiplicities.

###### Acknowledgements.

We thank the Formation and Evolution of Planetary Systems group for their
consistent support, Dr. Max Moe and Dr. Kimberly Ward-Duong for their
insightful advice, and an anonymous referee whose helpful reviews
significantly improved the manuscript.

## References

  * Bowler et al. (2015) Bowler, B. P., Liu, M. C., Shkolnik, E. L., & Tamura, M. 2015, ApJS, 216, 7
  * Burgasser et al. (2007) Burgasser, A. J., Reid, I. N., Siegler, N., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 427
  * Cortés-Contreras et al. (2017) Cortés-Contreras, M., Béjar, V. J. S., Caballero, J. A., et al. 2017, A&A, 597, A47
  * De Furio et al. (2019) De Furio, M., Reiter, M., Meyer, M. R., et al. 2019, ApJ, 886, 95
  * De Rosa et al. (2014) De Rosa, R. J., Patience, J., Wilson, P. A., et al. 2014, MNRAS, 437, 1216
  * Delfosse et al. (1999) Delfosse, X., Forveille, T., Beuzit, J. L., et al. 1999, A&A, 344, 897
  * Delfosse et al. (1998) Delfosse, X., Forveille, T., Perrier, C., & Mayor, M. 1998, A&A, 331, 581
  * Dieterich et al. (2012) Dieterich, S. B., Henry, T. J., Golimowski, D. A., Krist, J. E., & Tanner, A. M. 2012, AJ, 144, 64
  * Duchêne & Kraus (2013) Duchêne, G. & Kraus, A. 2013, ARA&A, 51, 269
  * Dupuy & Liu (2011) Dupuy, T. J. & Liu, M. C. 2011, ApJ, 733, 122
  * El-Badry et al. (2019) El-Badry, K., Rix, H.-W., Tian, H., Duchêne, G., & Moe, M. 2019, MNRAS, 489, 5822
  * Elliott et al. (2015) Elliott, P., Huélamo, N., Bouy, H., et al. 2015, A&A, 580, A88
  * Fischer & Marcy (1992) Fischer, D. A. & Marcy, G. W. 1992, ApJ, 396, 178
  * Fontanive et al. (2018) Fontanive, C., Biller, B., Bonavita, M., & Allers, K. 2018, MNRAS, 479, 2702
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Henry et al. (2006) Henry, T. J., Jao, W.-C., Subasavage, J. P., et al. 2006, AJ, 132, 2360
  * Janson et al. (2012) Janson, M., Hormuth, F., Bergfors, C., et al. 2012, ApJ, 754, 44
  * Jódar et al. (2013) Jódar, E., Pérez-Garrido, A., Díaz-Sánchez, A., et al. 2013, MNRAS, 429, 859
  * Marcy & Benitz (1989) Marcy, G. W. & Benitz, K. J. 1989, ApJ, 344, 441
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
  * Raghavan et al. (2010) Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1
  * Reggiani & Meyer (2013) Reggiani, M. & Meyer, M. R. 2013, A&A, 553, A124
  * Shan et al. (2015) Shan, Y., Johnson, J. A., & Morton, T. D. 2015, ApJ, 813, 75
  * Shvartzvald et al. (2016) Shvartzvald, Y., Maoz, D., Udalski, A., et al. 2016, MNRAS, 457, 4089
  * Ward-Duong et al. (2015) Ward-Duong, K., Patience, J., De Rosa, R. J., et al. 2015, MNRAS, 449, 2618
  * Winters et al. (2021) Winters, J. G., Charbonneau, D., Henry, T. J., et al. 2021, AJ, 161, 63
  * Winters et al. (2019) Winters, J. G., Henry, T. J., Jao, W.-C., et al. 2019, AJ, 157, 216