diff --git a/doc/omeas_heavy_mesons.qmd b/doc/omeas_heavy_mesons.qmd index 91774a64b..b1ef91910 100644 --- a/doc/omeas_heavy_mesons.qmd +++ b/doc/omeas_heavy_mesons.qmd @@ -20,6 +20,10 @@ csl: acta-ecologica-sinica.csl fontsize: 16pt --- +### Preamble and notation +
+ Show Content + `tmLQCD`can compute the correlators for the heavy mesons $K$ and $D$. The twisted mass formulation of the heavy doublet $(s, c)$ is such that $K$ and $D$ mix in the spectral decomposition @baron2011computing. In fact, the Dirac operator is not diagonal in flavor, @@ -54,8 +58,13 @@ In the following we use the twisted basis $\chi$ for fermion fields: \end{align} +
+ ## Correlators for the $K$ and $D$ +
+ Show Content + The required correlators are given by the following expectation values on the interacting vacuum ${\bra{\Omega} \cdot \ket{\Omega} = \braket{\cdot}}$: @@ -139,8 +148,13 @@ where ${\sigma_3 = 3. The action of $\sigma_1$ swaps the up and down flavor components of a spinor. ::: +
+ ## Wick contractions +
+ Show Content + Using the above remarks, we can write: @@ -208,7 +222,7 @@ Our correlator becomes: Upon a careful calculation for all values $i,j = 0,1$ we find, equivalently (using spacetime translational symmetry): -\begin{equation} +\begin{equation} \label{eq:C.hihj.Gamma1Gamma2} \begin{split} \mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) &= @@ -219,23 +233,53 @@ Upon a careful calculation for all values $i,j = 0,1$ we find, equivalently (usi (S_u)^\dagger (x|0) \gamma_5 \Gamma_1 \right] -\\ -&= -\sum_{y,z} -\operatorname{Tr} -\left[ - (S_h)_{f_i f_j} (x+y|y) - \Gamma_2 \gamma_5 - (S_u)^\dagger (x+z|z) - \gamma_5 \Gamma_1 - \delta_{yz} -\right] -\, . \end{split} \end{equation} This is a generalized case of eq. (A9) of @PhysRevD.59.074503. +
+ + +## Stochastic approximation of the correlators + +### Index dilution + +
+ Show Content + + +We now make the following remark. If we want to invert numerically the system $D_{ij} \psi_{j} = \eta_{i}$ (where the $i,j$ indices include all internal indices), we have: + +\begin{equation} +D_{ij} \psi_{j} = \eta_{i} \, \implies +\psi_{i} = S_{ij} \eta_{j} \, . +\end{equation} + +If $\eta_i \eta_j^{*} = \delta_{ij}$ we have: + +\begin{equation} +S_{ij} = \psi_{i} \eta_{j}^{*} \, . +\end{equation} + +--- + +We can also use **index dilution** in order to select the components we want. In fact, if we define: $\eta_i^{(a)} = \eta \delta_{i}^{a}$, with $\eta^{*} \eta =1$ we have: + +\begin{equation} +S_{ab} = S_{ij} \delta_{i}^{a} \delta_{j}^{b} = +S_{ij} (\eta^{*} \delta_{i}^{a}) (\eta \delta_{j}^{b}) = +\psi_{i}^{(b)} (\eta_{i}^{a})^{*} += \eta^{(a)} \cdot \psi^{(b)} \, . +\end{equation} + +
+ +### Stochastic expression of the correlators + +
+ Show Content + We now approximate the propagator using stochastic sources. Additionally, we use: @@ -246,7 +290,7 @@ Additionally, we use: \begin{equation} \eta^{(\alpha)}_{\beta, c} = \eta_c \, \delta^\alpha_\beta \,\, , \, \, - \eta^\dagger_c \eta_c = 1 \, . + \eta_c \otimes \eta^\dagger_c = \mathbb{1}_{N_c \times N_c} \, . \end{equation} - Flavor dilution: sources have an additional flavor index $\phi$, such that their flavor component different from the index vanish. The value of the flavor component however is the same, it changes only its position in the doublet: @@ -267,99 +311,63 @@ Additionally, we use: \end{equation} +Therefore, we can approximate the correlators of eq. \eqref{eq:C.hihj.Gamma1Gamma2} above with: -Therefore, we can use spin dilutions to rephrase the correlator in a form which will turn out to be convenient later -($c$ is the color index): - -\begin{equation} -\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) -= -[(S_h)_{f_i f_j}]_{\alpha_1 \beta_1} (x|0) -[\eta^{(\alpha_2)}_{\beta_1}]_c -(\Gamma_2 \gamma_5)_{\alpha_2 \alpha_3} -[{(\eta^\dagger)}^{(\alpha_3)}_{\beta_2}]_c -[(S_u)^\dagger]_{\beta_2 \alpha_4} (x|0) -(\gamma_5 \Gamma_1)_{\alpha_4 \alpha_1} -\, . -\end{equation} - - -We now define our spinor propagators. -If $\eta^{(\beta, \phi)}$ is the diluted source: - -\begin{align} -& (D_{\ell/h})_{\alpha_1 \alpha_2} (x|y) ({\psi}_{\ell/h}^{(\beta, \phi)})_{\alpha_2} (y) -= (\eta^{(\beta, \phi)})_{\alpha_1} (x) -\\ -& \, \implies \, -(\psi_{\ell/h}^{(\beta, \phi)})_{\alpha_1} (x) -= -(S_{\ell/h})_{\alpha_2 \alpha_1} (x | y) -\eta^{(\beta, \phi)}_{\alpha_2} (y) -= -(S_{\ell/h})_{\alpha_2 \alpha_1} (x | 0) -\eta^{(\beta, \phi)}_{\alpha_2} (0) -\\ -& \, \implies \, -(\psi_{\ell/h}^{(\beta, \phi)})^{*}_{\alpha_1} (x) -= -(\eta^{(\beta, \phi)})^{*}_{\alpha_2} (y) -(S_{\ell/h}^\dagger)_{\alpha_2 \alpha_1} (x | y) -= -(\eta^{(\beta, \phi)})^\dagger_{\alpha_2} (0) -(S_{\ell/h}^\dagger)_{\alpha_2 \alpha_1} (x | 0) -\end{align} - -This means that for our matrix of correlators we have to do ${4_D \times 2_f \times 2_{h,\ell}} = 16$ inversions. - -::: {.remark} -1. For the light doublet `tmLQCD` computes only $S_u$, which is obtained with $(\psi_h^{(\beta, f_0)})_{f_0}$. This is the only propagator we need. -2. For the heavy propagator, we can access the $(i,j)$ component of $S_h$ with $(\psi_h^{(\beta, f_j)})_{f_i}$. -::: - -Our correlator is given by the following expectation value -(no summation on flavor indices): \begin{equation} \begin{split} \mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) &= -\langle -[(\psi_h^{(\alpha_2, f_j)})_{f_i}]_{\alpha_1}(x) +\left( + \eta^{(f_i, \alpha_1)}(0) + \cdot + \psi_h^{(f_j, \alpha_2)}(x) +\right) (\Gamma_2 \gamma_5)_{\alpha_2 \alpha_3} -[(\psi_\ell^{(\alpha_3, f_0)})_{f_0}^\dagger]_{\alpha_4} (x) +\left( + \eta^{(0, \alpha_3)}(0) + \cdot + \psi_u^{(0, \alpha_4)}(x) +\right)^{*} (\gamma_5 \Gamma_1)_{\alpha_4 \alpha_1} -\rangle -\\ -&= -\braket{ -(\psi_\ell^{(\alpha_3, f_0)})_{f_0}^\dagger (x) -\cdot -(\gamma_5 \Gamma_1) -\cdot -(\psi_h^{(\alpha_2, f_j)})_{f_i}(x) -} -\, -(\Gamma_2 \gamma_5)_{\alpha_2 \alpha_3} \\ -&= -\mathcal{R}^{\alpha_3 \alpha_2} (\Gamma_2 \gamma_5)_{\alpha_2 \alpha_3} +&= +(\psi_h)_{f_i, \alpha_1}^{(f_j, \alpha_2)}(x) +(\Gamma_2 \gamma_5)_{\alpha_2 \alpha_3} +(\psi_u)_{\alpha_3}^{(0, \alpha_4)}(x)^{*} +(\gamma_5 \Gamma_1)_{\alpha_4 \alpha_1} \end{split} \end{equation} -More explicitly: - +In the end we have: + \begin{equation} -\begin{split} -\Gamma_2=1 &\implies \mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) -= \mathcal{R}^{00}+\mathcal{R}^{11}-\mathcal{R}^{22}-\mathcal{R}^{33} -\\ -\Gamma_2=\gamma_5 &\implies \mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) -= \mathcal{R}^{00}+\mathcal{R}^{11}+\mathcal{R}^{22}+\mathcal{R}^{33} -\end{split} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) += +[\Gamma_2 \gamma_5 (\psi_u)^{(0, \alpha_2)}(x)^{*}]_{\alpha_1} +[\gamma_5 \Gamma_1 (\psi_h)_{f_i}^{(f_j, \alpha_1)}(x) ]_{\alpha_2} \end{equation} - + +In our case $\Gamma_{1,2} = 1,\gamma_5$. Since we use the chiral basis ${\Gamma_{1,2}^{*} = (\Gamma_{1,2}^{T})^\dagger = \Gamma_{1,2}}$. +Therefore we can equivalently write the correlator as: + +\begin{equation} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) += +[\Gamma_2 \gamma_5 (\psi_u)^{(0, \alpha_2)}(x)]^{*}_{\alpha_1} +\cdot +[\gamma_5 \Gamma_1 (\psi_h)_{f_i}^{(f_j, \alpha_1)}(x) ]_{\alpha_2} \, +\end{equation} + +where $\cdot$ is the dot product in color space (NOTE: it complex-conjugates the 1st vector). + + + +
+ + + diff --git a/meas/correlators.c b/meas/correlators.c index 9bfe25db5..245e42020 100644 --- a/meas/correlators.c +++ b/meas/correlators.c @@ -40,7 +40,6 @@ #include "gettime.h" - #define TM_OMEAS_FILENAME_LENGTH 100 /****************************************************** @@ -253,6 +252,15 @@ void light_correlators_measurement(const int traj, const int id, const int ieo) return; } + +void spinor_dirac_array(su3_vector* phi, spinor psi){ + phi[0] = psi.s0; + phi[1] = psi.s1; + phi[2] = psi.s2; + phi[3] = psi.s3; +} + + /****************************************************** * * This routine computes the correlator matrix (), @@ -292,7 +300,7 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, #endif /* heavy-light correlators variables */ - // sign change bilinear^\dagger, with Gamma = 1,gamma_ + // sign change bilinear^\dagger, with Gamma = 1,gamma_5 double eta_Gamma[2] = {1.0, -1.0}; // even-odd spinor fields for the light and heavy doublet correlators @@ -473,9 +481,10 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, for (size_t beta = 0; beta < 4; beta++) { // spin dilution index for (size_t F = 0; F < 2; F++) { // flavor dilution index for (size_t i_f = 0; i_f < 2; i_f++) { // flavor index of the doublet - eo_source_spinor_field_spin_diluted_oet_ts(arr_eo_spinor[0][db][beta][F][0][i_f], - arr_eo_spinor[0][db][beta][F][1][i_f], t0, - beta, sample, traj, seed_i); + eo_source_spinor_field_spin_diluted_oet_ts( + arr_eo_spinor[0][db][beta][F][0][i_f], + arr_eo_spinor[0][db][beta][F][1][i_f], t0, + beta, sample, traj, seed_i); } } } @@ -498,15 +507,15 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, // stored temporarely in the propagator spinors (used as dummy) if (F==0){ mul_one_pm_itau2_and_div_by_sqrt2( - arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], - arr_eo_spinor[0][1][beta][F][i_eo][0], phi1, +1.0, - VOLUME / 2); + arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], + arr_eo_spinor[0][1][beta][F][i_eo][0], phi1, +1.0, + VOLUME / 2); } if (F==1){ mul_one_pm_itau2_and_div_by_sqrt2( - arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], - phi2, arr_eo_spinor[0][1][beta][F][i_eo][1], +1.0, - VOLUME / 2); + arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], + phi2, arr_eo_spinor[0][1][beta][F][i_eo][1], +1.0, + VOLUME / 2); } for (size_t i_f = 0; i_f < 2; i_f++) { @@ -598,9 +607,9 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, // (c,s) --> [(1-i*tau_2)/sqrt(2)] * (c,s) // stored temporarely in phi1, phi2 mul_one_pm_itau2_and_div_by_sqrt2(phi1, phi2, - arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], -1.0, - VOLUME / 2); - + arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], -1.0, + VOLUME / 2); + assign(arr_eo_spinor[1][1][beta][F][i_eo][0], phi1, VOLUME / 2); assign(arr_eo_spinor[1][1][beta][F][i_eo][1], phi2, VOLUME / 2); } @@ -616,19 +625,19 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, } // now we switch from even-odd representation to standard - //for (size_t i_sp = 0; i_sp < 2; i_sp++) { // source or sink - for (size_t db = 0; db < 2; db++) { // doublet: light of heavy - for (size_t beta = 0; beta < 4; beta++) { // spin dilution index - for (size_t F = 0; F < 2; F++) { // flavor projection index - for (size_t i_f = 0; i_f < 2; i_f++) { // flavor index - convert_eo_to_lexic(arr_spinor[1][db][beta][F][i_f], - arr_eo_spinor[1][db][beta][F][0][i_f], - arr_eo_spinor[1][db][beta][F][1][i_f]); + for (size_t i_sp = 0; i_sp < 2; i_sp++) { // source or sink + for (size_t db = 0; db < 2; db++) { // doublet: light of heavy + for (size_t beta = 0; beta < 4; beta++) { // spin dilution index + for (size_t F = 0; F < 2; F++) { // flavor projection index + for (size_t i_f = 0; i_f < 2; i_f++) { // flavor index + convert_eo_to_lexic(arr_spinor[i_sp][db][beta][F][i_f], + arr_eo_spinor[i_sp][db][beta][F][0][i_f], + arr_eo_spinor[i_sp][db][beta][F][1][i_f]); + } } } } } - // } //MPI_Barrier(MPI_COMM_WORLD); if (g_proc_id == 0){ @@ -637,7 +646,7 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, //MPI_Barrier(MPI_COMM_WORLD); if (g_proc_id == 0){ - printf("Checkpoint 11\n"); + printf("Checkpoint 11.0\n"); } /* @@ -645,18 +654,21 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, https://arxiv.org/pdf/1005.2042.pdf) I can build the correlators of eq. (20) */ const int f0 = 0; // flavor index of the up + /* now we sum only over local space for every t */ + const int j_ts = t0-g_proc_coords[0]*T; // checkerboard index of the time at the source for (t = 0; t < T; t++) { - j = g_ipt[t][0][0][0]; + j = g_ipt[t][0][0][0]; // source and sink separated by "t" // dummy variables res = 0.0; respa = 0.0; resp4 = 0.0; - for (i = j; i < j + LX * LY * LZ; i++) { + + // light correlators for (size_t beta = 0; beta < 4; beta++) { // spin dilution spinor psi_u = arr_spinor[1][0][beta][f0][f0][i]; @@ -665,25 +677,47 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, respa += _spinor_prod_re(psi_u, phi); _gamma5(phi, phi); resp4 += _spinor_prod_im(psi_u, phi); - - for (size_t hi = 0; hi < 2; hi++) { - for (size_t hj = 0; hj < 2; hj++) { - for (size_t g1 = 0; g1 < 2; g1++) { - for (size_t g2 = 0; g2 < 2; g2++) { - double dum = 0.0; // dummy variable - - // heavy doublet spinor propagator - phi = arr_spinor[1][1][beta][hj][hi][i]; - if (g1 == 0) { // Gamma_1 = \mathbb{1} - _gamma5(phi, phi); + } + + // heavy correlators + for (size_t hi = 0; hi < 2; hi++) { + for (size_t hj = 0; hj < 2; hj++) { + for (size_t g1 = 0; g1 < 2; g1++) { + for (size_t g2 = 0; g2 < 2; g2++) { + double complex dum_tot = 0.0; + for (int alpha_1 = 0; alpha_1 < 4; alpha_1++){ + spinor psi_h = arr_spinor[1][1][alpha_1][hj][hi][i]; + if (g1 == 0){ // Gamma_1 = Id + _gamma5(psi_h, psi_h); } - _spinor_scalar_prod(dum, psi_u, phi); - - double sign_beta = +1.0; - if (g2 == 0 /* Gamma_2 == \mathbb{1} */ && beta >= 2) { - sign_beta = -1.0; + su3_vector psi_h_su3[4]; + spinor_dirac_array(&psi_h_su3[0], psi_h); + for (int alpha_2 = 0; alpha_2 < 4; alpha_2++){ + spinor psi_u_star = arr_spinor[1][0][alpha_2][f0][f0][i]; + if (g2 == 0){ // Gamma_2 = Id. NOTE: works because Gamma_2=Gamma_2* for Gamma_2=1,gamma_5 + _gamma5(psi_u_star, psi_u_star); + } + su3_vector psi_u_star_su3[4]; + spinor_dirac_array(&psi_u_star_su3[0], psi_u_star); + complex double dum_12 = 0.0; + _colorvec_scalar_prod(dum_12, psi_u_star_su3[alpha_1], psi_h_su3[alpha_2]); + dum_tot += dum_12; } - res_hihj_g1g2[hi][hj][g1][g2] += sign_beta * dum; + } + // if (g_proc_id == 0){ + // printf("dum_tot = %d %d %d %d %d %e %e \n", t, hi, hj, g1, g2, creal(dum_tot), cimag(dum_tot)); + // } + if (hi != hj){ + dum_tot *= -1; + } + if (g1==g2){ + res_hihj_g1g2[hi][hj][g1][g2] = creal(dum_tot); // correlators is real + } + else{ + if (g2==1){ + dum_tot *= -1; + } + res_hihj_g1g2[hi][hj][g1][g2] = cimag(dum_tot); // correlator is imaginary } } } @@ -808,7 +842,8 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, for (size_t hj = 0; hj < 2; hj++) { for (size_t g1 = 0; g1 < 2; g1++) { for (size_t g2 = 0; g2 < 2; g2++) { - fprintf(ofs, "%u %u %u %u 0 %e %e\n", hi, hj, g1, g2, C_hihj_g1g2[hi][hj][g1][g1][t0], 0.0); + printf("check matrix: %e\n", C_hihj_g1g2[hi][hj][g1][g2][3]); + fprintf(ofs, "%u %u %u %u 0 %e %e\n", hi, hj, g1, g2, C_hihj_g1g2[hi][hj][g1][g2][t0], 0.0); for (t = 1; t < g_nproc_t * T / 2; t++) { tt = (t0 + t) % (g_nproc_t * T); fprintf(ofs, "%u %u %u %u %d %e ", hi, hj, g1, g2, t, C_hihj_g1g2[hi][hj][g1][g2][tt]); @@ -904,11 +939,12 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, } void correlators_measurement(const int traj, const int id, const int ieo) { - //light_correlators_measurement(traj, id, ieo); // ??? maybe add a double check? does i1 --> light and i2 --> heavy? if (measurement_list[id].measure_heavy_mesons == 1) { const int i1 = 0, i2 = 1; heavy_correlators_measurement(traj, id, ieo, i1, i2); } + + light_correlators_measurement(traj, id, ieo); } diff --git a/offline_measurement.c b/offline_measurement.c index 95ea094e8..4911ecd14 100644 --- a/offline_measurement.c +++ b/offline_measurement.c @@ -293,9 +293,7 @@ int main(int argc, char *argv[]) if (g_proc_id == 0) { fprintf(stdout, "#\n# Beginning offline measurement.\n"); } - //printf("Ciao Simone4!\n"); meas->measurefunc(nstore, imeas, even_odd_flag); - printf("Ciao Simone5!\n"); } nstore += Nsave; } diff --git a/solver/dirac_operator_eigenvectors.h b/solver/dirac_operator_eigenvectors.h index b4ef212e1..6b7454533 100644 --- a/solver/dirac_operator_eigenvectors.h +++ b/solver/dirac_operator_eigenvectors.h @@ -265,6 +265,14 @@ int cyclicDiff(int a,int b, int period); conj((r).s3.c2) * (s).s3.c2; + +#define _colorvec_scalar_prod(proj,r,s)\ + (proj) = conj((r).c0) * (s).c0 + \ + conj((r).c1) * (s).c1; + \ + conj((r).c2) * (s).c2; + + + #define PROJECTSPLIT(p_plus,up_plus,col_proj,phi_o,phi_plus,col_phi)\ p_plus = 0; \ p_plus += conj(up_plus->s0.col_proj) * (phi_o->s0.col_phi); \