looking for proof or partial proof of determinant conjecture Math people:
I am looking for a proof of a conjecture I made.  I need to give two definitions. For distinct real numbers $x_1, x_2, \ldots, x_k$, define $\sigma(x_1, x_2, \ldots, x_k) =1$ if $(x_1, x_2, \ldots, x_k)$ is an even permutation of an increasing sequence, and $\sigma(x_1, x_2, \ldots, x_k) =-1$ if $(x_1, x_2, \ldots, x_k)$ is an odd permutation of an increasing sequence. For example, $\sigma(2, 1, 10, 8) = 1$ because $(2,1,10,8)$ is an even permutation of $(1,2,8,10)$, and $\sigma(2, 1, 8, 10) = -1$ because $(2,1,8,10)$ is an odd permutation of $(1,2,8,10)$.  For real $B$, $n\geq 1$ and distinct real numbers $\mu_1, \mu_2, \ldots, \mu_n, \gamma_1, \gamma_2, \ldots, \gamma_n$, let $M(B;\mu_1,\mu_2,\ldots,\mu_n;\gamma_1,\gamma_2,\ldots,\gamma_n)$ be the $n$-by-$n$ matrix defined by 
$$ M(B;\mu_1,\mu_2,\ldots,\mu_n;\gamma_1,\gamma_2,\ldots,\gamma_n)_{i,j}=\frac{\exp(-B\gamma_j)}{\mu_i+\gamma_j}+\frac{\exp(B\gamma_j)}{\mu_i-\gamma_j}.        $$
My conjecture is the following: if $n \geq 1$, $B \geq 0$, and $\mu_1, \mu_2, \ldots, \mu_n, \gamma_1, \gamma_2, \ldots, \gamma_n$ are distinct positive numbers with 
$0<\mu_1 < \mu_2 < \cdots < \mu_n$ and $0<\gamma_1 < \gamma_2 < \cdots < \gamma_n$ , then
$$\operatorname{sgn}(\operatorname{det}(M(B;\mu_1,\mu_2,\ldots,\mu_n;\gamma_1,\gamma_2,\ldots,\gamma_n))) = (-1)^{\frac{n(n+1)}{2}}
           \sigma(\mu_1, \mu_2, \ldots, \mu_n, \gamma_1, \gamma_2, \ldots, \gamma_n). $$
Of course $\operatorname{sgn}(x)$ is the sign of $x$, which is $1$, $-1$, or $0$. I have proven this is true for $n=1$ and $n=2$.  For $n$ between $3$ and $20$, I have run thousands of experiments in Matlab using randomly generated $\mu$'s and $\gamma$'s.  In a set of one thousand experiments, the conjectured equation will typically hold every single time, or might fail once or twice, with the determinant (with the wrong sign) being extremely small, so perhaps roundoff error is the culprit.
UPDATE: let $d(B)$ be the determinant of the matrix, where the other parameters should be clear from context.  $d(B)$ is an analytic function of $B$.  It suffices to show that $\frac{\partial^m d}{\partial B^m}$ has the desired sign at $B=0$ for all $m \geq 0$.  Unfortunately, the determinant of $\frac{\partial M}{\partial B}$ is not the same thing as $\frac{\partial^m d}{\partial B^m}$ (if it were, properties of Cauchy matrices would yield the desired conclusion).  Since the conjecture is true for  $n=1$, that means that the displayed formula for $M_{i,j}$ above, and all its derivatives with respect to $B$, have the same sign as $\mu_i - \gamma_j$ at $B=0$, and $M_{i,j}$ has that sign for all positive $B$.  I proved the conjecture for $n=2$, by computing the determinant of $M$ and its derivatives with respect to $B$ at $B=0$, and looking at the six possible orderings of $\mu_1, \mu_2, \gamma_1$, and $\gamma_2$ given the restrictions $\mu_1 < \mu_2$ and $\gamma_1 < \gamma_2$.  I had some help from Maple multiplying out, simplifying and factoring algebraic expressions.  I am trying to prove the general case by induction on $n$, expanding the determinant along the last row or column, but the determinants of the $n-1$-by-$n-1$ minors don't seem to necessarily have the ``right'' signs.
Thanks to some comments provided below, unless I am confused, the conjecture can be proven for $B=0$ and large positive $B$, for any $n$, using properties of Cauchy matrices.
 A: I have a proof if all the $\mu$'s are larger than all of the $\gamma$'s. I've spent a while trying to figure out how modify this proof to work for other orderings, but I'm giving up. This argument is inspired by a computation of the Caucy determinant by Doron Zielberger (skip to the 25 minute mark). 
We start with the identity
$$e^{\mu B} \int_B^{\infty} e^{- \mu t} (e^{\gamma t} + e^{-\gamma t}) dt =   \frac{e^{-B \gamma}}{\mu+\gamma} + \frac{e^{B \gamma}}{\mu-\gamma}.$$
Note that we need $\mu > \gamma$ in order for the integral to converge; that is why I can't extend this proof to any other ordering of the $\mu$'s and $\gamma$'s
So the determinant is
$$\det \left( e^{\mu_i B} \int_{t_i=B}^{\infty} e^{-\mu_i t_i} (e^{\gamma_j t_i} + e^{- \gamma_j t_i}) d t_i \right).$$
Note that we are using $n$ different integration variables -- one for each row of the matrix. By the multilinearity of the determinant, this is
$$\exp\left(\sum B \mu_i \right) \times \int_{t_1=B}^{\infty} \int_{t_2=B}^{\infty} \cdots \int_{t_n=B}^{\infty} \exp\left( - \sum \mu_i t_i \right) \det \left( e^{\gamma_j t_i} + e^{- \gamma_j t_i} \right) dt_1 dt_2 \cdots dt_n.$$
Note that permuting $(t_1, \ldots, t_n)$ only changes $ \det \left( e^{\gamma_j t_i} + e^{- \gamma_j t_i} \right)$ by a sign, but changes $\exp(\mu_i t_i)$ to $\exp(\mu_{\sigma(i)} t_i)$ for some permutation $\sigma$. Lumping together all $n!$ reorderings of the $t_i$, the integral is
$$\int_{B \leq t_1 \leq \cdots \leq t_n} \sum_{\sigma \in S_n} (-1)^{\sigma} \exp(\sum t_i \mu_{\sigma(i)}) \det \left( e^{\gamma_j t_i} + e^{- \gamma_j t_i} \right) dt_1 \cdots dt_n $$
$$= \int_{B \leq t_1 \leq \cdots \leq t_n} \det(e^{\mu_j t_i}) \det \left( e^{\gamma_j t_i} + e^{- \gamma_j t_i} \right) dt_1 \cdots dt_n.$$
We claim that the integrand is positive for $0 < t_1 < \cdots < t_n$ so the integral is positive, as desired. So we are done once we prove:
Lemma If $\mu_1 < \mu_2 < \cdots < \mu_n$ and $\gamma_1 < \gamma_2 < \cdots < \gamma_n$ and $t_1 < t_2 < \cdots < t_n$, then
$$\det \left( e^{\mu_j t_i} \right) \ \mbox{and} \  \det \left( e^{\gamma_j t_i} + e^{- \gamma_j t_i} \right)$$
are positive.
Proof It is enough to show that these determinants don't vanish anywhere in this range, since they then must have constant sign and it is easy to check that the correct sign is positive. If the first determinant vanished, then there would be some nonzero function $\sum a_j z^{\mu_j}$ which vanished at $z = e^{t_1}$, $e^{t_2}$, ... $e^{t_n}$. But then this is a "polynomial with real exponents" having $n$ nonzero terms, and $n$ postive roots, contradicting Descartes rule of signs. (I wrote out a proof of Descartes' rule of signs for real exponents here.)
Similarly, if the second determinant vanishes, then we get a real-exponent-polynomial $\sum b_j (z^{\gamma_j} + z^{-\gamma_j})$ with $2n$ terms and roots at the $2n$ points $e^{\pm t_i}$. Again, a contradiction.
A: Reports on some experiments: If all the $\gamma$'s and $\mu$'s are integers, then the expression is a Laurent polynomial in $e^B$ with integer coefficients. I decided to write $z$ for $e^B$ and test this case. Some basic Mathematica code for anyone who wants it
mm[x_, y_] := z^(-y)/(x + y) + z^y/(x - y)

dd[xlist_, ylist_] := Expand[Det[Table[mm[x, y], {x, xlist}, {y, ylist}]]]

So dd[{1,3,5},{2,4,6}] will return the determinant for $\mu= \{ 1,3,5 \}$ and $\gamma = \{ 2,4,6 \}$.
DoIt[xlist_, ylist_] :=
  If[Length[Union[xlist, ylist]] < Length[xlist] + Length[ylist], 
    Null,
    With[{r = NSolve[{z > 1, dd[xlist, ylist] == 0}, z]}, 
      If[Length[r] > 0, 
        Print[{xlist, ylist, r}], 
        Null]]]

DoItVerbose[xlist_, ylist_] :=
  If[Length[Union[xlist, ylist]] < Length[xlist] + Length[ylist], 
    Print["Collision"],
    With[{r = NSolve[{z > 1, dd[xlist, ylist] == 0}, z]}, 
      If[Length[r] > 0, 
        Print[{xlist, ylist, r}], 
        Print[{xlist, ylist, "No roots"}]]]

DoItLots[l_, t_] := 
 Do[DoIt[Table[Random[Integer, {1,100}], {l}], 
   Table[Random[Integer, {1,100}], {l}]], {t}]

In other words, DoIt[xlist, ylist] searches for a counterexample with xlist and ylist for $\mu$ and $\gamma$, printing the counterexample if it finds it and otherwise doing nothing. DoItVerbose is similar but it will tell you the manner in which it failed. DoItLots[l,t] will try $t$ times to find a counterexample, using $l \times l$ determinants.
I tried $500$ times to find a counterexample with $3 \times 3$ determinants, and $20$ times each for $4 \times 4$, $5 \times 5$, ..., $10 \times 10$. Around here, the birthday paradox starts being an issue, meaning that we waste a lot of time with $\gamma$ and $\mu$ having repeated values, so I would do something more intelligent before going further.
I also decided to test the suggestion in the comments above that the determinant might have positive coefficients as a power series in $B$. 
 DoItPowerSeries[xlist_, ylist_] :=
   If[Length[Union[xlist, ylist]] < Length[xlist] + Length[ylist], 
    Null,
    If[Length[
        Union[Sign[
         CoefficientList[
           Series[dd[xlist, ylist] /. z -> E^b, {b, 0, 10}],
           b]]]] > 1, 
      Print[{xlist, ylist}], 
      Null]]

DoItPSLots[l_, t_] := 
  Do[DoItPowerSeries[Table[Random[Integer, {1, 1000}], {l}], 
    Table[Random[Integer, {1, 1000}], {l}]], {t}]

In other words, expand the power series out to $b^{11}$ and check whether all the coefficients have the same sign. (By Igor Rivin's argument, we know what the sign will be for $b$ sufficiently large, so it was easier to test for "all same sign" then to explicitly do the $\pm 1$ computation.) 100 tests each for $3 \times 3$ and $4 \times 4$ without a counterexample.  I also ran the same tests with clearing out the denominator (so the matrix entries are now $1/(\mu-\gamma) + z^{2 \gamma}/(\mu+\gamma)$ and substituting $z=w+1$; again, all coefficients appear to have the same sign.
A: I don't have anything like a complete answer. BUT, consider the case where $B$ is very large in absolute value (and assume the $\gamma$s are all positive or all negative). Then (assuming that $\gamma$s are positive, and $B \gg 0$) your matrix elements are approximately 
$$\frac{\exp(B\gamma_j)}{\mu_i - \gamma_j}.$$ The exponentials factor out of the determinant to give you a multiplicative factor of $\exp(B \sum \gamma_j),$ and what remains is a Cauchy matrix of which the determinant can be evaluated explicitly, and it seems like your conjecture is correct. Similarly, when $B\ll 0$ you get $\exp(-B\sum\gamma_j)$ times a (different) Cauchy determinant. You should check that your conjecture is correct then too. If, for some reason, you knew that your matrix was always nonsingular as $B$ varies you would be golden, but this is not clear. It is clear that the determinant is zero whenever some $\gamma_i$ is equal to some $\gamma_j,$ and also clear that the determinant has a pole whenever some $\mu_i$ equals some $\mu_j,$ so if the determinant were a rational function, I think you would know that these were all the zeros and poles, but it is not. Still, this might be close.
NOT true for negative $B$ For $2\times 2$ matrix with $\mu_1 = 1, \mu_2=3, \gamma_1 = 2, \gamma_2=4$ the determinant is positive for positive $B$ and very negative $B$ but is negative in a range of negative $B$ (but not VERY negative, which makes me think that the conjecture is actually false, and the countexemples the OP found might be real counterexamles. In fact, this is also borne out by a (sort-of) closed form for the determinant: 
closed form Do the dumbest thing possible and use the multilinearity of the determinant, to get $$ D = \sum_{\mbox{subsets of $1, \dots, n$}} \exp(B\sum \mathbf{\gamma}_S) C_{\mathbf{\mu}, \mathbf{\gamma}_S},$$
where $\mathbf{\gamma}_S$ is the set of $\gamma$s with the ones in $S$ negated, while $C_{\mathbf{\mu}, \mathbf{\gamma}_S}$ is the Cauchy determinant, with $\mu$s and the some-signed-flipped $\gamma$s.
