Updating Bounds on $R$-Parity Violating Supersymmetry from Meson Oscillation Data

We update the bounds on $R$-parity violating supersymmetry originating from meson oscillations in the $B^0_{d/s}$ and $K^0$ systems. To this end, we explicitly calculate all corresponding contributions from $R$-parity violating operators at the one-loop level, thereby completing and correcting existing calculations. We apply our results to the derivation of bounds on $R$-parity violating couplings, based on up-to-date experimental measurements. In addition, we consider the possibility of cancellations among flavor-changing contributions of various origins, e.g. from multiple $R$-parity violating couplings or $R$-parity conserving soft terms. Destructive interferences among new-physics contributions could then open phenomenologically allowed regions, for values of the parameters that are naively excluded when the parameters are varied individually.


Introduction
Several years of operation of the LHC have (as yet) failed to reveal any conclusive evidence for physics beyond the Standard Model (SM) [1]. On the contrary, experimental searches keep placing ever stronger limits on hypothesized strongly [2-5] and even weakly-interacting [6] particles in the electroweak-TeV range. While this situation tends to leave the simpler models in an uncomfortable position, for the so-called "CMSSM" see for example Ref. [7], it also advocates for a deeper study of more complicated scenarios, satisfying the central motivations of the original paradigm but also requiring more elaborate experimental investigations for testing.
Softly-broken supersymmetric (SUSY) extensions of the SM [8,9] have long been regarded as a leading class of candidates for the resolution of the hierarchy problem [10], as well as a possible framework in view of understanding the nature of dark matter or the unification of gauge-couplings. The simplest of such models, the Minimal Supersymmetric Standard Model (MSSM), has thus been the focus of numerous studies in the past decades. An implicit ingredient of the usual MSSM is R-parity (R p ) [11], a discrete symmetry related to baryon and lepton number. In addition to the preservation of these quantum numbers, R p is also invoked in order to justify the stability of the lightest SUSY particle, leaving it in a position of a dark-matter candidate [12].
Despite its attractive features, R p conservation is not essential to the phenomenological viability of a SUSY model. R p violation (RpV) -see [13,14] for reviews -is viable as well; simply a different discrete (or gauge) symmetry is required [15][16][17][18]. It also leads to a distinctive phenomenology which is relevant to LHC searches [19,20].
With experimental constraints now coming from both low-energy physics and the highenergy frontier, it seems justified to give the RpV-phenomenology a closer look, beyond the tree-level or single-coupling approximations that are frequently employed in the literature.
In this paper, we consider the most general RpV-model with minimal superfield content. The superpotential of the R p -conserving MSSM is thus extended by the following terms [21]: where Q,Ū ,D, L,Ē denote the usual quark and lepton superfields, · is the SU (2) L invariant product and ε abc is the 3-dimensional Levi-Civita symbol. The indices i, j, k refer to the three generations of flavor, while a, b, c correspond to the color index. We note that symmetryconditions may be imposed on the parameters λ ijk and λ ijk without loss of generality: λ ijk = −λ jik , λ ijk = −λ ikj . The first three sets of terms of Eq.(1.1) violate lepton-number and the last set of terms violates baryon-number.
The superpotential of Eq.(1.1) contains several sources of flavor-violation, in both the lepton and the quark sectors. Such effects are steadily searched for in experiments, placing severe bounds on the parameter space of the model. The impact of lepton-flavor violating observables on the RpV-MSSM has been discussed extensively in the literature, see e.g. . In the quark sector, observables such as leptonic B-decays or radiative b → s transitions [47][48][49] have been considered. Here, we wish to focus on neutral-meson mixing observables, ∆M K , ∆M d , ∆M s , for K 0 , B 0 d and B 0 s mesons, respectively. Such observables have been discussed in the R-parity conserving [50,51] as well as in an RpV context in the past [47,[52][53][54][55][56][57][58][59]. Yet, diagrams beyond the tree-level and box contributions as well as sfermion or RpV-induced mixings have been routinely ignored. The purpose of this paper consists in addressing these deficiencies and proposing a full one-loop analysis of the meson-mixing observables in the RpV-MSSM. These values are in excellent agreement with the SM computations [61,62], resulting in tight constraints on new physics contributions.
For the K 0 −K 0 system, the Particle Data Group [63] combines the experimental measurements as: ∆M exp K = (0.5293 ± 0.0009) · 10 −2 ps −1 . (1.3) Despite the precision of this result, constraints from K 0 −K 0 mixing on high-energy contributions are considerably relaxed by the large theoretical uncertainties due to long-distance effects. Historically, estimates of the latter have been performed using the techniques of large N QCD -see e.g. Ref. [64] -while lattice QCD collaborations such as [65] are now considering the possibility of evaluating these effects in realistic kinematical configurations.
Ref. [66] settles for a long-distance contribution at the level of (20 ± 10)% of the experimental value, and we follow this estimate below. Concerning short-distance contributions, Ref. [67] performed a NNLO study of the charm-quark loops, resulting in a SM estimate of ∆M SM, Short Dist. K = (0.47 ± 0.18) · 10 −2 ps −1 .
The computation of the meson oscillation parameters is usually performed in a lowenergy effective field theory (EFT), where short-distance effects intervene via the Wilson coefficients of dimension 6 flavor-changing (∆F = 2) operators [68]. This procedure ensures a resummation of large logarithms via the application of the renormalization group equations (RGE) from the matching high-energy (e.g. electroweak) scale down to the low-energy (meson-mass) scale where hadronic matrix elements should be computed [69]. In this work, we calculate the contributions to the Wilson coefficients arising in the RpV-MSSM up to oneloop order. The λ couplings of Eq.(1.1) already generate a tree-level diagram. Going beyond this, at one-loop order, diagrams contributing to the meson mixings involve both R-parity conserving and R-parity violating couplings. These are furthermore intertwined via RpVmixing effects stemming for example from the bilinear term µ i H u · L i . Our analysis goes beyond the approximations that are frequently encountered in the literature. We also find occasional differences with published results, which we point out accordingly.
In the following section, we= present the general ingredients of our full one-loop analytical calculation of the Wilson coefficients of the ∆F = 2 EFT (effective field theory) in the RpV-MSSM, referring to the appendices where the exact expressions are provided. In Section 3, we discuss our implementation of these results employing the public tools SPheno [70,71], SARAH [72][73][74][75][76][77], FlavorKit [78] and Flavio [79]. Finally, numerical limits on the RpV-couplings are presented in a few simple scenarios, before a short conclusion.

Matching conditions for the ∆F = EFT of the RpV-MSSM
We consider the ∆F = 2 EFT relevant for the mixing of (d i d j )-(d j d i ) mesons -d i corresponds to the down-type quark of ith generation (d, s or b). The EFT Lagrangian is written as where we employ the following basis of dimension 6 operators: The superscripts (a, b = 1, 2, 3) refer to the color indices when the sum is not trivially contracted within the fermion product. We have employed the usual four-component spinor notations above, with P L,R denoting the left-and right-handed projectors.
The Wilson coefficients C i ,C i associated with the operators of Eq.(2.5) in the Lagrangian of the EFT -Eq.(2.4) -are obtained at high-energy by matching the d idj → d jdi amplitudes in the EFT and in the full RpV-MSSM. We restrict ourselves to the leading-order coefficients (in a QCD/QED expansion) on the EFT-side. On the side of the RpV-MSSM, we consider only short-distance effects, i.e. we discard QCD or QED loops. Indeed, the photon and gluon are active fields in the EFT, so that a proper processing of the corresponding effects would require a NLO matching procedure. Furthermore, both tree-level and one-loop contributions are considered in the RpV-MSSM: we stress that this does not induce a problem in power-counting, as the tree-level contribution is a strict RpV-effect, so that R p -conserving (or violating) one-loop amplitudes are not (all) of higher QED order. Numerically speaking, one possibility is that the tree-level is dominant in the Wilson coefficients, in which case, the presence of the one-loop corrections does not matter. This case is essentially excluded if we consider the experimental limits on the meson-oscillation parameters. If, on the contrary, the tree-level contribution is of comparable (or subdominant) magnitude with the one-loop amplitudes, then the electroweak power-counting is still satisfied. Yet, one-loop contributions that are aligned with the tree-level always remain subdominant.
For our calculations in the RpV-MSSM, we employ the Feynman 't Hooft gauge [80] and dimensional regularization [81,82]. For reasons of consistency with the tools that we employ for the numerical implementation, DR-renormalization conditions will be applied. However, in the results that we collect in the Appendix, the counterterms are kept in a generic form, which allows for other choices of renormalization scheme. We apply the conventions where the sneutrino fields do not take vacuum expectation values. 7 Moreover, the λ couplings of Eq.(1.1) are defined in the basis of down-type mass-states, i.e. a CKM matrix appears when the second index of λ connects with an up-type field, but not when it connects to a downtype field [52]. Mixing among fields are considered to their full extent, including left/right and flavor squark mixings, charged-Higgs/slepton mixing, neutral-Higgs/sneutrino mixing, chargino/lepton mixing and neutralino/neutrino mixing. The details of our notation and the Feynman rules employed can be found in Appendix A. As a crosscheck, we performed the calculation using two different approaches for the fermions: the usual four-component spinor description and the two-component description [85].
On the side of the EFT, the operators of Eq.(2.5) each contribute four tree-level Feynman diagrams to the d idj → d jdi amplitude. Half of these contributions are obtained from the other two by an exchange of the particles in the initial and final states: as the dimension 6 operators are symmetrical over the simultaneous exchange of both d i 's and both d j 's, we may simply consider two diagrams and double the amplitude. The two remaining diagrams correspond to an (s ↔ t)-channel exchange. We exploit these considerations to reduce the number of diagrams that we consider on the side of the RpV-MSSM to only one of the s/t-channels.
The tree-level contribution to the d idj → d jdi amplitudes is due to the λ couplings of 7 For the general rotation to this basis see Ref. [83]. See also Ref. [84] for a discussion of this in terms of physics at the unification scale. Eq.(1.1). It involves a sneutrino exchange where, however, sneutrino-flavor and sneutrino-Higgs mixing could occur. The appearance of RpV contributions at tree-level complicates somewhat a full one-loop analysis: one-loop contributions indeed depend on the renormalization of the d idj -sneutrino vertex (and of its external legs). In principle, one could define this vertex 'on-shell', i.e. impose that one-loop corrections vanish for on-shell d i , d j external legs -while the counterterm for the sneutrino field is set at momentum p 2 = M 2 K,B 0. In such a case, one could restrict oneself to calculating the box-diagram contributions to d idj → d jdi . However, in any other renormalization scheme, self-energy and vertex-correction diagrams should be considered. Yet, if the λ couplings contributing at tree-level are small, the impact of the vertex and self-energy corrections is expected to be limited, since these contributions retain a (at least) linear dependence on the tree-level λ . These contributions are symbolically depicted in Fig.1.
One-loop diagrams contributing to d idj → d jdi include SM-like contributions (box diagrams with internal u, c, t quarks, W and Goldstone bosons), 2-Higgs-doublet-model-like contributions (box diagrams with internal u, c, t quarks, charged-Higgs bosons and possibly W or Goldstone bosons), R p -conserving SUSY contributions (box diagrams with chargino/scalarup, neutralino/sdown or gluino/sdown particles in the loop) and RpV-contributions (selfenergy and vertex corrections, box diagrams with sneutrino/quark, slepton/quark, lepton/squarks, neutrino/squark or quark/squark internal lines). The RpV-driven mixing further intertwines these contributions, so that the distinction among e.g. the R p -conserving chargino/scalar-up and RpV lepton/scalar-up boxes becomes largely superfluous. For all these contributions, with exception of the self-energy diagrams on the external legs, we neglect the external momentum, as it controls effects of order m d i,j , which are subdominant when compared to the  momentum-independent pieces of order M W or M SUSY . Yet, when a SM-fermion f appears in the loop, some pieces that are momentum-independent still come with a suppression of order m f /M W,SUSY . We keep such pieces even though they could be discarded in view of the previous argument.
The diagrams of Fig.1 are calculated in Appendix B (tree-level contribution), Appendix C (d i -quark self-energies), Appendix D (scalar self-energy) and Appendix E (vertex corrections). While we go beyond the usual assumptions employed to study the ∆F = 2 Wilson coefficients in the RpV-MSSM, it is possible to compare the outcome of our calculation to partial results available in the literature. First, in the limit of vanishing RpV-parameters, we recover the well-known results in the R p -conserving MSSM, which are summarized in e.g. the appendix of Ref. [50]. Then, RpV-contributions from the tree-level and box-diagram topologies have been presented in Ref. [47] in the no-mixing approximation. Taking this limit and neglecting further terms that are not considered by this reference, we checked that our results coincided, with the exception of the coefficient c λ LR of Ref. [47] (a piece of the contribution to C 5 ). Transcripted to our notations, the result of Ref. [47] reads: while we obtain: The mismatch lies in the prefactor and the sfermion chiralities. Another class of λ boxes involving an electroweak charged current has been considered in the no-mixing limit in Ref. [54]. There, we find agreement with our results. As self-energy and vertex corrections have not been considered before, the opportunities for comparison are more limited. Still, we checked that the scalar self-energies were consistent with the results of Ref. [86]. Finally, our results can be controlled in another fashion, using the automatically generated results of public tools: we detail this in the following section.

Numerical implementation and tools
In order to determine limits from the meson oscillation measurements on the parameter space of the RpV-MSSM, we establish a numerical tool implementing the one-loop contributions to the ∆F = 2 Wilson coefficients and deriving the corresponding theoretical predictions for ∆M K,d,s . To this end, we make use of the Mathematica package SARAH [72][73][74][75][76][77] to produce a customized spectrum generator based on SPheno [70,71,87]. SPheno calculates the complete supersymmetric particle spectrum at the one-loop order and includes all important two-loop corrections to the neutral scalar masses [88].
The routines performing the calculation of flavor observables are generated through the link to FlavorKit [78]. FlavorKit makes use of FeynArts/FormCalc [89][90][91] to calculate the leading diagrams to quark and lepton flavor violating observables. For the meson mass differences, the tree-level and box diagrams as well as the double-penguin contributions are included per default. However, as parameters within SPheno are defined in the DR scheme, it is in principle necessary to implement the self-energy and vertex corrections. We added the vertex corrections via PreSARAH [78], which enables the implementation of new operators into FlavorKit within certain limits. As the scalar self-energies cannot be generated in this fashion, we incorporated these by hand.
The Wilson coefficients computed by FlavorKit and PreSARAH at the electroweak matching scale are stored in analytical form in the Fortran output of FlavorKit. We compared these expressions with our results of the previous section; we found explicit agreement in almost all cases -and adapted the code to match our results in the few cases where it proved necessary. 8 After the Wilson coefficients at the electroweak matching scale are computed, further steps are necessary in order to relate them to the observables ∆M K,d,s . The FlavorKit output includes a theoretical prediction for these observables, however the hadronic input parameters are more up-to-date in the more recently-developed code Flavio [79], which shares an interface with FlavorKit using the FLHA standards [92]. We hence use Flavio to process the Wilson coefficients as calculated by FlavorKit. First, the Wilson coefficients must be run to a low-energy scale using the QCD RGE's of the EFT [69]. In the case of the K 0 −K 0 system, the impact of the charm loop is sizable [67]: we upgraded the NLO coefficient η cc coded within Flavio to the NNLO value 1.87 (76) [67] and η ct = 0.496(47) [93]. For consistency, the charm mass in the loop functions is set to the MS value m c (m c ) 1.28 GeV. Then, the hadronic dynamics encoded in the dimension 6 operators must be interpreted at low-energy in the form of hadronic mixing elements: this step gives rise to "bag-parameters", which are evaluated in lattice QCD. Here, Flavio employs the bag parameters of Ref. [94] for the K 0 −K 0 system and of Ref. [95] for the B 0 d −B 0 d and B 0 s −B 0 s systems. In addition, the CKM matrix elements within Flavio are derived from the four inputs |V us |, |V ub |, |V cb | and γ. We set these to the fit-results of Ref. [63]: |V us | 0.22506, |V ub | 3.485 · 10 −3 , |V cb | 4.108 · 10 −2 and γ 1.236. Moreover, we changed the B 0 d decay constant to a numerical value of 186 MeV [96]. Finally, we added the observable ∆M K to Flavio (based on pre-included material) and made sure that the predicted SM short-distance prediction was consistent with the theoretical SM estimate given by Ref. [67].
A quantitative comparison of the predicted ∆M K,d,s with the experimental results of Eqs.(1.2) and (1.3) requires an estimate of the theoretical uncertainties. The Wilson coefficients have been obtained at leading order, which implies higher-order corrections of QCD-size. In the case of the SM-contributions, large QCD logarithms are resummed in the evolution of the RGEs between the matching electroweak scale and the low-energy scale. 8 In rare cases, we identified seemingly minor -but numerically important -differences between our computation and the FlavorKit code, namely in a few tree-level contributions to C5 (which should be absent), as well as inC2,3 and C2,3 for a few one-loop box diagrams. We fixed those appearances in the code as well as the relative sign between tree and one-loop contributions after correspondence and cross-checking with the However, for the new-physics contributions, further logarithms between the new-physics and the electroweak scale could intervene -FlavorKit computes the new-physics contributions to the Wilson coefficients at the electroweak scale, hence missing such logarithms. Therefore, the higher-order uncertainty is larger for contributions beyond the SM and can be loosely es- , where µ N P and µ EW represent the new-physics and electroweak scales, respectively. Further sources of uncertainty are the RGE evolution in the EFT and the evaluation of hadronic matrix elements. For the SM matrix elements, the uncertainties on η cc , η ct and η tt are of order 30% [67], 10% [93] and 1% [97], respectively, leading to a large SM uncertainty in ∆M K and a smaller one in ∆M d,s . For the K 0 −K 0 system, the bag-parameters are known with a precision of ∼ 3% in the case of B K | for the short-distance contribution to ∆M K . As explained above, we will employ the estimate of Ref. [

Numerical results
We are now in a position to study the limits on RpV-parameters that are set by the mesonoscillation parameters. However, it makes limited sense to scan blindly over the RpV-MSSM parameter space imposing only constraints from the ∆M 's. Comparable analyses of all the relevant observables for which experimental data is available would be necessary. We will thus restrict ourselves to a discussion of the bounds over a restricted number of parameters and in a few scenarios. The input parameters that we mention below correspond to the SPheno input defined at the M Z scale.
We first consider the case where no explicit source of flavor violation appears in the R pconserving parameters. The flavor transition is thus strictly associated to the CKM matrix or to the RpV-effects. The latter can intervene in several fashions: • Flavor violation in the λ couplings could also intervene at the loop-level only. This happens when, for instance, one product of the form λ mnI λ * mnJ or λ mIn λ * mJn is non-zero -again, (I, J) corresponds to the valence quarks of the meson; m and n are internal to the loop.
• Finally, the flavor transition can be conveyed by the λ couplings, in which case it appears only at the loop level in the ∆M 's. Possible coupling combinations include λ m12 λ m23 , λ m12 λ m13 or λ m13 λ m23 .
Below, we first consider these three cases separately, before we investigate possible interferences between tree-and loop-level generated diagrams for several non-zero λ couplings. However, we avoid considering simultaneously non-zero LQD andŪDD couplings: then, discrete symmetries no longer protect the proton from decay, so that the phenomenology would rapidly come into conflict with associated bounds. Still, we note that some diagrams contributing to the meson mixing parameters would combine both types of couplings: these are also provided in the appendix.
Then, flavor transitions can also be mediated by R p -conserving effects. In this case, flavor violation could originate either in the CKM matrix, as in the Minimal Flavor Violation scenario [98], or in new-physics parameters, such as the soft squark bilinear and trilinear terms. We briefly discuss possible interferences with RpV-contributions.
For simplicity, we consider only the case of real λ ( ) and disregard the bilinear R-parity violating terms (though they are included in our analytical results in the appendix).

Bounds on a pair of simultaneously non-zero LQD couplings
Scenario   Figure 3: Constraints from the ∆M 's on scenarios with RpV-mediated flavor violation contributing at tree-level, as a function of the sneutrino mass. The plots on the left correspond to the upper limit on positive λ · λ ; those on the right to lower limits on negative λ · λ combinations. The green, orange, red and purple colors represent regions within [0, 1σ], [1σ, 2σ], [2σ, 3σ] and > 3σ bounds, respectively. The experimental central value is exactly recovered on the black lines.

Tree Level Contributions
Let us begin with the case where only two LQD couplings are simultaneously non-vanishing and contribute to the ∆M 's at tree-level. For doing so, we choose a spectrum of the form of an effective SM at low mass, where we have fixed the squark, higgsino and gaugino masses to 2 TeV, while varying all the slepton masses simultaneously in the range 0.2 − 2 TeV. The important parameter values are listed in the first line of Table 1. In addition, the stop trilinear coupling A t , of order 3 TeV, is adjusted so that the lighter Higgs mass satisfies m h ≈ 125 GeV. We also considered several other scenarios, listed in Table 1, e.g. involving lighter charged Higgs or lighter squarks of the third generation, but the general properties of the constraints remained qualitatively unchanged. This indicates that new-physics R pconserving contributions (mediated by the CKM matrix) to the meson oscillation parameters are relatively small, in view of the uncertainties. All the input is defined at the electroweak scale, so that we can discuss the various classes of RpV-contributions to the ∆M 's without the blurring effect due to the propagation of flavor-violation via RGE's between a high-energy scale and the electroweak scale.
In Fig. 3, we present the limits set by ∆M d , ∆M s and ∆M K on the tree-level flavor violating contributions. The plots in the first column are obtained for a positive product λ · λ , while those in the second column correspond to negative λ · λ . For each observable, the most relevant λ · λ combination, leading to a tree-level contribution, was selected. The individual sub-figures depict the extension of the 0, 1, 2, 3 σ regions in the plane defined by the corresponding flavor-violating λ · λ product and the slepton mass. The colors in Fig. 3 are chosen such that purple regions are excluded at three standard deviations or more; red regions are excluded at ≥ 2 σ -which is the limit that we apply later on, in order to decide whether a point in parameter space is excluded or allowed experimentally; the orange regions correspond to a prediction of the ∆M within 1 and 2 σ; finally, the green areas are consistent with the experimental measurement within 1 σ, while the black curves reproduce the central values exactly. Experimental and theoretical uncertainties are added in quadrature to define the total uncertainty U tot = U 2 theo + U 2 exp . In the case of ∆M K , the theoretical uncertainties from long-distance and short-distance contributions are also combined quadratically. Since experimentally one cannot tell apart the two mass eigenstates of B 0 d/s , we simply consider the absolute value of ∆M d/s in our evaluation. When we plot ∆M d,s , this feature may result in a doubling of the solutions for the central value or of the 1 σ-allowed regions, such as in the upper-left and middle-left plots of Fig. 3. For K 0 , instead, the mass ordering, and hence the sign of ∆M K is known.
The limits that we obtain on the λ couplings contributing at tree-level are relatively tight. In the scenarios of Fig.3, the 2σ bounds read approximately: where we assume that only one lepton flavor, namely i, has non-vanishing RpV-couplingstherefore the bounds only depend on the mass of the corresponding sneutrinoν i . Alternatively, with degenerate sneutrinos, we could sum over the index i on the left-hand side of Eq. (4.8). Limits on these products of couplings have been presented in Ref. [99] for a SUSY mass of 100 GeV and in [59] for a mass of 500 GeV -as explained above, our limits can be confronted to the bounds applying on i λ i13 λ i31 , etc., in these references. In comparison, the bounds that we obtain in Fig.3 are somewhat stronger, at least by a factor ∼ 3. This result should be put mainly in the perspective of the reduction of the experimental uncertainty in the recent years.

1-Loop Contributions to Flavor Transition
Next, we turn to the case where a pair of LQD couplings mediate the flavor transition only at the loop-level and we focus on coupling combinations of the form λ mnI λ * mnJ or λ mIn λ * mJn (with I, J the valence quarks of the meson). In principle we could consider other combinations, such as λ mnI λ * mnJ , λ mnI λ * mñJ , λ mIn λ * mJn or λ mIn λ mJñ (with m =m, n =ñ). However, either the associated contributions are CKM suppressed or they would require several λ · λ products to be simultaneously non-zero or non-degenerate scalar / pseudoscalar sneutrino fields. We thus restrict ourselves to the two types mentioned above. For these, we note that the limits are independent of the flavor m of the slepton field. The spectrum is described in the third row of Table 1. In this context, RpV-effects in ∆M 's are dominated by diagrams involving the comparatively light (charged or neutral) sleptons. We thus concentrate on these below. We can distinguish two types of contributions: • If one of the pair of non-vanishing LQD couplings is one of those involved for the treelevel exchange diagram -i.e. if it contains the two flavor indices of the valence quarks of the meson -we find that quark self-energy corrections on the tree-level diagram can be comparable to or even dominant over box contributions.
• If neither of the non-vanishing LQD couplings participates in the tree-level diagrams, box diagrams are the main contributions.
This difference impacts both the magnitude of the resulting bounds and their dependence on the slepton mass, as we shall see below. In Fig.4, we consider non-vanishing λ 121 λ 123 , λ 112 λ 113 and, finally, λ 113 λ 123 . In these cases, the box diagrams dominate over the fermionic self-energy corrections. For each scenario, the limits from the ∆M 's essentially originate in one of the three observables ∆M d , ∆M s or ∆M K . The corresponding limits approximately read: where ml i denotes the mass of the degenerate sneutrinos and charged sleptons. Here, we note that the mass dependence of the form (λ · λ ) 2 < c · m 2 differs from that appearing when the RpV-contribution intervenes at tree-level. It is characteristic of the leading RpV-diagrams in the considered setup, corresponding to the box formed out of two charged sleptons and two up-type quarks in the internal lines and to the box consisting of two sneutrinos and two down-type quarks: these diagrams roughly scale as (λ ·λ ) 2 /m 2 . As a consequence, the limits for positive and negative λ · λ products are comparable. In addition, the bounds on λ · λ now scale about linearly with the sparticle mass.
Expectedly, the limits are much weaker in these box-dominated scenarios than in the case where the flavor transition appears at tree-level. Refs. [54,55,59] presented limits on the corresponding coupling-combinations for a sfermion mass of 100 or 500 GeV. The bounds that we derive are of the same order. Corresponding scenarios are displayed in Fig.5, where ∆M B d , ∆M Bs and ∆M K are plotted against λ 131 · λ 133 , λ 132 · λ 133 and λ 121 · λ 122 , respectively. The bounds have a comparable scaling to that appearing in the scenario with tree-level sneutrino exchange, but the constraints are far weaker. At 2 σ: Due to the inclusion of the missing and obviously relevant self-energy diagrams, the bounds that we report are accordingly tighter than in the literature [54,55,59].
In Table 2, we compile the 2 σ bounds on λ ·λ products that we derive for 1 TeV sleptons in the scenario SUSY-RpV(a) of Table 1 (the limits depend only weakly on the chosen scenario). Table 2: Compilation of the latest bounds on relevant couplings of LQD operators, coming from the considered meson oscillation observables. These limits were established with the spectrum defined in the row SUSY-RpV(a) of Table 1, with slepton and sneutrino masses of 1 TeV. The precise 2 σ boundary obviously depends on the sign of the non-vanishing λ · λ product: we always apply the most conservative (weakest) limit. In the list of couplings, the comment "(T)/(S)/(B)" indicates that the coupling product is dominated by a treelevel/quark self-energy/box contribution. "N/A" means that we did not identify upper-limits on the couplings below 4π (a rough limit from perturbativity considerations). Above the horizontal line, the non-vanishing coupling combinations select right-handed external quarks. Below this line, the external quarks are left-handed. The scaling with the sneutrino/slepton mass is roughly quadratic for all λ · λ products that contain both valence flavors in (at least) one of the non-vanishing λ , linear otherwise: see more precise explanation in the main body of the text. Some combinations contribute to two observables, such as λ i13 λ i32 , relevant for both ∆M d and ∆M s . In such a case, the square brackets identify the weaker limit. In this list, the pairs λ · λ are taken non-zero only one at a time and, in particular, for a unique (s)lepton flavor i. As explained above, the scaling with the slepton/sneutrino mass depends on the choice of non-vanishing λ : essentially quadratic if at least one of the non-vanishing λ contains both valence-flavors of the decaying meson, linear otherwise. One of the ∆M 's is usually more sensitive to a specific λ · λ product than the other two. etc. We proceed with our analysis and now consider baryonic RpV, i.e. non-zeroŪDD couplings. The corresponding RpV-effects appear only at the radiative level and are dominated by box diagrams. Contrarily to existing analyses [47], we always consider heavy gluinos (as indicated by the current status of LHC searches), so that the associated diagrams generally remain subdominant. In this setup, three classes of diagrams compete: (1) boxes including two squarks and two quarks in internal lines, which scale like (λ · λ ) 2 , (2) boxes including two quarks, one squark and a W -boson, which scale like λ ·λ but involve a CKM-suppression and a quark-chirality flip, and (3) similarly boxes with two squarks, one quark and a chargino, which scale like λ ·λ . The matter of the chirality flip can be easily understood as only righthanded quarks couple via λ but only left-handed quarks couple to a W . Therefore, such diagrams with an internal W line are mostly relevant when the internal quark line involves a top-quark. As to the boxes with an internal chargino line, we also find that such contributions are mainly relevant for an internal stop line: indeed, the higgsino contribution scales with the Yukawa coupling, hence is suppressed for squarks of first or second generation. In addition, the gaugino contribution relies on left-right squark mixing, which we keep negligible for squarks of the first and second generation -making the assumption that the trilinear soft terms are proportional to the Yukawa couplings [8].

Bounds on a pair of simultaneously non-zeroŪDD couplings
In Fig. 6, we present the 1, 2 and 3 σ limits on coupling combinations allowing for internal (s)top lines. The relevant right-handed squarks are assumed to be mass-degenerate. The regime with small λ couplings is dominated by the box diagrams involving W bosons and top quarks in the internal lines. We find that, for low mass values, this contribution scales with the squark mass in an intermediate fashion between linear and quadratic, because of the finite top mass effects. These effects largely vanish for squark masses above O(1 TeV) and we then recover the scaling with λ ·λ m 2 q . The supersymmetrized version of the W boxes, i.e.
boxes with internal charginos, are also contributing with a scaling of λ ·λ m 2 q . However their impact w.r.t. the W boxes is always reduced. At large values of the couplings and for light squarks, the purelyŪDD-mediated diagrams appear to be the most relevant, scaling with (λ ·λ ) 2 m 2 q -in analogy to the slepton box-diagrams with non-vanishing LQD coupling -so that the bounds on λ · λ show a roughly linear dependence with the squark mass. Then, for both large |λ · λ | and heavier quarks, the W -mediated diagrams and these purelyŪDD boxes can be of comparable magnitude, hence lead to interference structures. This interplay between various contributions brings about a non-trivial mass dependence of the bounds on the λ couplings, with both constructive as well as destructive effects between the individual amplitudes. The plots for negative λ · λ couplings perfectly illustrate this fact, in particular in the case of ∆M s . Beyond this interference regime, at sufficiently large squark masses, the contribution from the UDD box with an internal W-line eventually supersedes the pure UDD amplitude. Since the bounds on the individual coupling combinations do not scale with a simple power law in mq R , we refrain from showing approximate expressions as we did in the scenarios with flavor-violation of LQD-type.
In Fig. 7, by contrast, the choice of non-vanishing λ couplings does not allow for internal (s)top lines. Thus the RpV-diagrams with mixed W /squark or chargino/quark internal lines are suppressed, and the scaling of the limits from meson-oscillation parameters is closer to linear. In addition, the 2 σ bounds are somewhat milder than in the previous case and roughly symmetrical for positive and negative λ · λ products. Thus, in this case, we extract the approximate bounds onŪ 1DiDj coupling pairs: Given that the scaling of the bounds on λ · λ pairs decidedly depends on the specific choice of couplings, we refrain from showing a compilation table as Table 2 for the LQD couplings, since it would only be representative of a specific SUSY spectrum.

Competition among LQD-driven contributions
Bounds on individual RpV-coupling products may be misleading, in the sense that several RpV-effects could cancel one another. In fact, the decomposition along the line of the low-energy flavors provides likely-undue attention to these specific directions of RpV, while the latter have no deep specificity from the high-energy perspective. In particular, RGE's are expected to mix the various flavor-directions of non-vanishing RpV-couplings, while the boundary condition at, say, the GUT scale, has no particular reason for alignment with the low-energy flavor directions [84,100]. Obviously however, the relevant directions in flavor space are highly model-dependent and we have no particular suggestion to make from the low-energy perspective of this work. Instead, we simply wish to illustrate the possibility of allowed directions with large RpV-couplings. To this end, we allow for two non-vanishing λ · λ coupling products and investigate the limits originating in the ∆M measurements.
If we consider Figs. 3 and 5, the tree-level diagram for λ i31 · λ i13 = O(10 −6 ) and the RpV-box for λ i31 · λ i33 = O(10 −4 ) -implying a hierarchy λ i13 /λ i33 = O(10 −2 ) -naively contribute to ∆M d by amplitudes of comparable magnitude. Whether these contributions can interfere destructively clearly depends on the form of the amplitudes but also on the sign of the non-vanishing couplings. In Fig. 8, we complete the results from Figs. 3 and 5 by now allowing for three non-vanishing couplings. In practice, we set the slepton/sneutrino mass to 1 TeV and keep one LQD coupling to a constant value: λ 131 = 0.01, λ 132 = 0.1, or λ 121 = 0.1. Then, we vary two independent λ , our choice depending again on the valence quarks of the considered ∆M . However, we stress that this procedure in fact opens three non-trivial λ · λ directions, so that the game is somewhat more complex than just playing one contribution  Table 1, with slepton/sneutrinos of 1 TeV. As in the previous plots, the color code reflects the level of tension between our predictions and the experimental measurements.
versus the other.
As expected, in the plots of Fig. 8, the interplay of various RpV-contributions opens funnel-shaped allowed regions for comparatively large values of the LQD couplings, highlighting the possibility of destructive interferences. We note that, considering that the treelevel and radiative contributions do not necessarily have the same scaling with respect to the slepton/sneutrino mass, the 'allowed angle' depends on the sfermion spectrum. Of course, the choice of parameters falling within the allowed funnels appears to be fine-tuned from the perspective of this work, but might be justified from a high-energy approach. On the other hand, constructive interferences lead to the 'rounded edges' observed in some of the plots.
As mentioned earlier, we will not consider the interplay of LQD-andŪDD-couplings, since such scenarios are of limited relevance without a quantitative analysis of the proton decay rate. On the other hand, our discussion in this subsection points to the relevance of considering a full evaluation of the ∆M 's (and other observables), when considering RpVscenarios beyond the simplistic one-coupling-dominance approach.

Competition between flavor violation in the R-parity conserving and R-parity violating sectors
RpV-couplings are not the only new sources of flavor violation in SUSY-inspired models.
In fact, the large number of possible flavor-violating parameters of the R p -conserving soft-SUSY-breaking Lagrangian is often perceived as a weakness for this class of model, known as the SUSY Flavor Problem. In particular, the soft quadratic mass-terms in the squark sector m 2 Q,Ū ,D and the trilinear soft terms A U,D are matrices in flavor-space that are not necessarily aligned with the flavor-structure of the Yukawa/CKM matrices. In this case, flavor-violation is generated in L−L, R−R (form 2 ) or L−R (for A) squark mixing. Correspondingly, flavorchanging-neutral gluinos or neutralinos, as well as new flavor-changing chargino couplings, could contribute to ∆M K,d,s in e.g. diagrams of the form of Fig. 2, (b-d) -see e.g. Ref. [55]. Here, we wish to illustrate the potential interplay of R p -conserving and RpV flavor violation. In particular, we note that the presence of flavor-violating effects in RpV-couplings would likely mediate flavor-violation in the squark sector via the RGE's [100].
We will focus on R p -conserving flavor-violation in the quadratic squark mass parameters m 2 ij , where we assume the diagonal terms to be degenerate for squarks of left-handed and right-handed type (for simplicity): m 2D = m 2 Q ≡ m 2 . Flavor-violation in the trilinear soft terms would lead to comparable effects at the level of the meson-oscillation parameters. However, large A-terms easily produce new (e.g. color-and charge-violating) minima in the scalar potential, that lead to instability of the usual vacuum, with possibly short-time tunnelling. In fact, we find that such stability considerations 9 typically constrain the A-terms much more  Table 1, with the slepton/sneutrino mass at 1.5 TeV. The flavor-violating quadratic soft mass parameters in the squark sector, m 2 ij , are chosen to be degenerate for left-handed and right-handed squarks. The color code follows the same conventions as before.  Table 1. The color code is unchanged compared to previous plots.
efficiently than the ∆M 's.
In Fig. 9, we allow for non-vanishing m 2 13 , m 2 23 or m 2 12 , simultaneously with non-zero λ 113 λ 131 , λ 123 λ 132 and λ 112 λ 121 . The former induce contributions to ∆M d , ∆M s and ∆M K through R p -conserving squark mixing, while the latter provide RpV tree-level contributions to the same ∆M 's. The parameters are set to the scenario SUSY-RpV(a) of Table 1, with the slepton/sneutrino mass at 1.5 TeV. In analogy with the results of section 4.3, we observe that R p -conserving and RpV contributions may interfere destructively or constructively. Thus, allowed funnels with comparatively large values of the RpV-couplings open. In particular, we note that a tiny m 2 12 is sufficient for relaxing limits from ∆M K , while the typical values of m 2 13 and m 2 23 affecting ∆M d and ∆M s are significantly larger. A similar analysis can be performed with RpV of theŪDD-type. This is shown in Fig. 10.
In this subsection, we have stressed that the limits originating from meson-oscillation parameters are quite sensitive to the possible existence of flavor-violating sources beyond that of the RpV-couplings. A full analysis of these effects thus appears necessary when testing a complete model.

Conclusions
In this paper, we have analyzed the meson-mixing parameters ∆M d,s and ∆M K at the full one-loop order in the RpV-MSSM. In particular, we have completed earlier calculations in the literature, in which only tree-level and box diagrams were usually considered. We also performed a numerical study based on our results and employing recent experimental and lattice data. The tighter limits that we derive -as compared to older works -illustrate the improvement of the precision in experimental measurements, but also the relevance of some of the new contributions that we consider. In particular, the interplay of SM-like and LQD-type flavor-violation modifies the scaling of the bounds with the sneutrino/slepton mass for a whole class of couplings. Finally, we have emphasized the possibility of interference effects amongst new sources of flavor violation, either exclusively in the RpV-sector or in association with R p -conserving squark mixing. While the appearance of allowed directions with comparatively large couplings largely intervenes as a fine-tuned curiosity in the low-energy perspective of our work, it also stresses the relevance of a detailed analysis of the observables when considering a complete high-energy model, since accidental relations among parameters could affect the picture of low-energy limits. [102]. A parameter point is deemed unstable on cosmological time-scales, and therefore ruled out, if the mean tunnelling time is smaller than 21.7% of the age of the Universe. Acknowledgements H. K. D. and Z.S.W. are supported by the Sino-German DFG grant SFB CRC 110 "Symmetries and the Emergence of Structure in QCD". M.E.K. is supported by the DFG Research Unit 2239 "New Physics at the LHC". The work of F. D. was supported in part by the MEINCOP (Spain) under contract FPA2016-78022-P, in part by the "Spanish Agencia Estatal de Investigación" (AEI) and the EU "Fondo Europeo de Desarrollo Regional" (FEDER) through the project FPA2016-78022-P, and in part by the AEI through the grant IFT Centro de Excelencia Severo Ochoa SEV-2016-0597. V.M.L. acknowledges support of the BMBF under project 05H18PDCA1. We thank the authors of FlavorKit, and in particular Avelino Vicente, for vivid exchanges and invaluable assistance in cross-checking our results with the FlavorKit routines. Likewise, we thank David Straub for very helpful correspondence concerning Flavio. We also thank Martin Hirsch and Toby Opferkuch for useful discussions. Z.S.W. thanks the IFIC for hospitality, and thanks the COST Action CA15108 for the financial support during his research stay at the IFIC.

A.1 Mixing matrices
• The squark mass matrices mix left-and right-handed components. We define the masseigenstates in terms of a unitary rotation of the gauge/flavor-eigenstates: Here, U α (resp. D α ) represents the scalar-up (resp. sdown) mass state with mass m Uα (resp. m Dα ). Summation over the generation index f is implicit.
• R-parity violation leads to a mixing of charged-Higgs and slepton fields. We define the mass-eigenstates H ± α with mass m Hα as: • Similarly, the neutral Higgs mass-states involve both the doublet-Higgs, H 0 ; in the CPviolating case, CP-even and CP-odd components mix as well.
S α denotes the mass-eigenstate associated with the mass m Sα .
• The charged winosw + ,w − , higgsinosh + u ,h − d and lepton fields e f L ,ē f R define the chargino sector. For the mass m χ ± k , the associated eignstate is given by: (A.15) • The violation of R-parity also mixes neutrino and neutralino states. The eigenstate with mass m χ 0 k reads:

A.2 Feynman rules
Here, we list the various couplings that intervene are relevant in our calculation.

A.3 Loop-functions
The loop functions relevant for our computations are • C 0 (p 1 , p 2 , m 1 , m 2 , m 3 ) = −16π 2 ı d D k Explicit expressions for these functions in the limit of vanishing external momenta can e.g. be found in Ref. [103].

B Tree level contributions
The tree-level contribution to the d idj → d jdi amplitudes corresponds to the topology of Fig.1a and is mediated by a sneutrino internal line. It generates the following terms in the EFT: where the couplings g

C d i − d j self-energy contributions
Loop corrections on the external d-fermion legs are determined by the LSZ reduction. Defining the matrix of renormalized d i − d j self energies as:Σ ij (p /) =Σ ij L (p /)P L +Σ ij R (p /)P R = P LΣ ij L (p /) + P RΣ ij R (p /), we derive the contribution to the EFT: We list below the contributions to the self-energies.

C.2 Vector/fermion loop
The vector/fermion pair (S/f ) is summed over the following list of particles: • W /up: Eq.(A.24).

C.3 Counterterm
Defining the generic d-mass counterterm δm d ji = δm L d ji P L + δm R d ji P R as well as the d-wavefunction counterterm δZ d ji = δZ L d ji P L + δZ R d ji P R , we arrive at the following contribution:

D Sneutrino-Higgs self-energies
We assume that the tadpoles (Higgs, gauge bosons) vanish, which supposes certain relations at the loop-level between vevs and tree-level parameters. Then, defining the renormalized neutral-scalar self-energy matrixΣ S αβ , we derive the following contribution to the EFT: (D.50) The various contributions to the neutral-scalar self-energies are listed below.

D.1 Scalar
This contribution is summed over the scalarS, taking value in the following list of particles: • scalar-ups: couplings from Eq.(A.41). 3 colors contributing.

D.2 Vector
The vector V belongs to the following list of particles: • W's: couplings from Eq.(A.39).

E Vertex corrections
The vertex corrections to the EFT are obtained as: receives the contributions listed below.

E.1 Scalar/fermion loop with cubic scalar coupling
List of particles for the scalar/fermion triplet (S k , S l /f ):

E.3 Vector/fermion loop with scalar-vector coupling
The vector/fermion triplet (V k , V l /f ) takes the following values:

E.4 Vector/fermion loop with scalar-fermion coupling
The vector/fermion triplet (V /f k , f l ) takes the following values: •

E.6 Counterterms
The counterterm contribution −ıV Sαd j d i [CT ] reads: is the counterterm to the Yukawa coupling and δλ R f ji = δλ L f ji * is the counterterm to the λ coupling.

F Box diagrams
Here, we collect the box-diagram contributions to the d idj → d jdi amplitude. The results are listed according to the topologies of Fig.2.
Case 2: f k colour-singlet List of particles: Case 3: f k colour-triplet