Hi Tom, thanks for this interesting reference. There are probably many more cases concerning  the impact of surfaces on enhancing the dynamics of surface versus bulk molecules. Years ago Jagadeesh Bharatram has observed such an effect when he vistited our lab (Chemical Physics Letters 404 (2005) 177–181). He found for lethicin covering the walls of nanopores, that the molecular dynamics in the surface controlled layers are less hindered compared to those in the bulk. It could be an interesting endeavor to consider all those phenomena in a unified context. But I believe, that for our manuscript, this would be beyond its scope. 

We certainly agree that we did not use the conventional expression for entropy in our simulations. This is why we state that we use crude models. In particular, we state beginning with line 197 "The free energy change Δ𝐴= Δ𝑈 − 𝑇 Δ𝑆 is determined from crude models of the internal energy change Δ𝑈 = 𝑭 Δ𝑹 defined by the repulsive net force 𝑭 exerted from all neighboring particles on the particle at stake and the length |Δ𝑹| of the jump to the jump to the next cell, the temperature 𝑇, and the entropy change Δ𝑆.“ Moreover, beginning with line 210 we write "It (the entropy) is estimated by the sum of the step lengths…“. We could have also approximated W in S = k_B lnW by the logarithm of the number of cells the particle is free to jump to. Instead, we used the sum of jump distances, which for our Moore neighborhood can be argued to approximate W (but not the logarithm) apart from some scaling factor. This we have chosen to approximate the configurational entropy for the discrete states in our simulation instead of the textbook formula S = -k_B S (P lnP). In fact, we used the sum of distances, because we are dealing with discrete configurations and the configurations on the square grid differ, so that S = k_B lnW does not strictly apply. Our choice may not be the best one, but our crude approximation exhibits the essential features of entropy: The distance sum is zero, if there is only one possible configuration, and it grows with the number of accessible configurations. For purpose of calculating jump probability this suffices. Now, one may argue, that the resulting jump probabilities are biased due to an incorrect substitution of S and that that bias leads to the asymmetry parameters we observed. For sure, the jump probabilities introduced in this way enhanced the observed asymmetry parameter, so that we detected it in the first place. But subsequently we found at much longer computation time when using a much higher number of steps that we could pick the destination cell in vacancy diffusion at random without calculating jump probabilities. This means that that principal existence of asymmetry in three-site exchange does not depend on the introduction of certain thermodynamic parameters or particle interactions in the simulation. So, the in the end, the way we approximated entropy is not critical to the conclusions reached in the manuscript, but it helped in the discovery process with the computer experiments on vacancy diffusion. That is why we describe our entropy substitute in the manuscript and, also, why we considered using gas diffusion calculations that are independent of our entropic substitution to verify the asymmetry parameter predictions.  Indeed, the exchange asymmetry observed for vacancy diffusion was also observed using the completely different gas-diffusion algorithm. This other program modelled elastic collisions between gas atoms starting with a Maxwell-Boltzmann distribution of initial velocities without considering the free energy.

Thank you for your persistence and for taking the time to educate us on computer modelling. Obviously, our expertise is not in this area. Nevertheless, we hoped to gain some insight into the symmetry properties of relaxation-exchange maps in NMR in this way.

R         Thank you for clarifying some of the questions. Regarding the technical issues, I still have questions.



R         1)  I am obviously not part of the NMR community, but I still think there has to be an equation for the asymmetry parameter which can be understood by everybody. Do I understand correctly that you average the quantity given in (5) over all jumps from 2->3 during the simulation, but you give different weights depending on the initial (or final, or both?) position of the particle (whether it is close to an active site or not?)

A         There are no different weights assigned to particles in different pools. The quantity Mj is a concentration. In the simulations it is the number of particles in pool j. The quantity k_ij denotes the rate of transitions from pool j to pool i. The program counts the number of particles passing from pool  j to i and assigns that to k_ij M_j. The denominator of (5) is the total number of jumps from one pool to another including jumps within one pool (k_ii Mi). So, the asymmetry parameter defined in (5) is the number difference between forward and backward jumps divided by the total number of jumps. The total number of jumps is the sum of all differential jumps, i. e. the sum over all k_ij M_j. It is calulated in the program after completion of each simulation run and has been verified by comparison with the initially specified number of jumps.

R         2) Thanks for clarifying the details of the off-lattice simulation. It is of course ok and common practice to use simulation units and not SI units. Your simulation units are apparently defined in terms of the mass m of the particles, the particle diameter σ (I assume it is one in your units), and the energy (I assume you set kT=1 when setting up the Maxwell-Boltzmann velocity distribution). This defines the time unit τ=(m s² /kT)^1/2  . In these units, your time step is Δt = 1 τ, which is very large. In molecular dynamics simulations, typical values are Δt = 10^-3 τ or less, and having a too large time step can have a severe impact on the results.

A         We agree that the time step has been too large and as a result has introduced a false image of asymmetry by an unintentional overlap of the hard circles and by possibly skipping collisions. The simulation has now been tested at shorter time steps to determine if such an error was introduced. We find that the asymmetry parameter gets smaller as the time step is decreased. Furthermore, we are looking into the reviewer’s suggestion of modifying the gas diffusion code to accomplish event driven calculations.

R         On the other hand, you do not really describe a Molecular Dynamics simulation, it is rather another type of (off-lattice) Monte Carlo simulation. For example, your collisions preserve the energy of the two colliding particles, but not their momentum. As a Monte-Carlo simulation, it does not preserve detailed balance, and therefore, again, it is not surprising that the results also break detailed balance.

A         The gas phase simulation is a common time-driven elastic hard circle model with walls rather than periodic boundaries. An initial distribution of speeds is generated and the particles are given a random initial direction of travel. Immediately after a collision, new velocities and deflection angles for the two particles are determined from conservation of momentum and kinetic energy by:

As mentioned in the previous response, the asymmetry parameter decreases with shorter time step suggesting that detailed balance is recovering for the shorter times.  Including event driven dynamics into this code, as also mentioned above, should help in the study of detailed balance violation.

R         3) Regarding the comment "We found nonzero asymmetry parameters also when choosing the destination cell for a jump from all vacant neighbor cells at random." I would like to note that this algorithm also breaks detailed balance. In order to maintain detailed balance, you have to choose a jump randomly from all neighbor cells and then reject the move if the neighbor cell is filled. Rejecting means the particle stays where it is and does not move at all, but the follow-up configuration still counts for the overall statistics.  But I suppose the authors are aware of this, since they also state " Choosing an algorithm that satisfies detailed balance by construction precludes testing detailed balance" (a statement to which I fully agree.)

A         Thank you for clarifying! We reran the calculations for Fig. 6e choosing randomly from all neighbor positions, free and occupied. Indeed, the asymmetry parameter produces noise more than one order of magnitude lower than the values observed with zero probability assigned to randomly chosen jumps to occupied positions. We now understand, that by introducing a probability to jumps, detailed balance is violated. So, we are driving the imbalance by our vacancy diffusion algorithm.

R         4) Regarding the comment on currents: Cyclic exchange, as long as it is persistent and does not average to zero on the long run, would also count as current in the  definition of thermal equilibrium.

A         ok.

R         5) Regarding thermalization: This not only means thermalization with respect to velocities. Asking whether a system thermalizes is the same as asking whether a system reaches thermal equilibrium in the above sense, i.e., a stationary state without stationary (macroscopic or microscopic) currents.  Testing this is generally difficult in a simulation. The velocity distribution usually approaches the Maxwell-distribution very quickly, but other quantities usually equilibrate much more slowly. In your simulations, you could for example test the system without force terms and check whether the system ever reaches the Boltzmann distribution with respect to positions. In the absence of any forces, the particles should be uniformly distributed in the pore. This is probably not the case.

A         What happens in thermal equilibrium is at the heart of the question. Detailed balance demands no currents. Experimental evidence (Song 2000) and tested theory (Brownstein 1977) indicate the existence of diffusion eigenmodes of fluids confined to pores, which describe spatially oscillating distributions of nuclear magnetization components in thermal equilibrium and which decay in distinct ways under the impact of diffusion and the boundary geometry. Considering, that three-site exchange probes Fick’s second law, we interpret our observed currents, to be a pore-resonance effect on translational motion which obeys the diffusion equation.

Conclusion

Thanks to the lessons learned from the expert advice of the reviewer of our manuscript we understand that our manuscript needs major revision. The main point is that the reported computer simulations do not indicate a violation of the principle of detailed balance. Therefore, we did not find any indication that three-site NMR exchange maps are asymmetric in thermodynamic equilibrium, and the observed asymmetry needs to be attributed to experimental or data processing (e.g. inverse Laplace transformation) deficiencies. A second outcome of the computer simulations remains without major revision. This is that the observed particle dynamics obey the diffusion equation and are linked to diffusion eigenmodes with the consequence that diffusion eigenmodes may possibly be driven by external stimuli like the violin bow enforcing resonance vibration of a Chladni plate. This is a technically interesting perspective as heterogeneous catalysis may possibly be enhanced be oscillating electrical, magnetic, or mechanical (ultrasound) fields. Pending the editor’s decision, we would revise our manuscript accordingly.

Thank you for your persistence and for taking the time to educate us on computer modelling. Obviously, our expertise is not in this area. Nevertheless, we hoped to gain some insight into the symmetry properties of relaxation-exchange maps in NMR in this way.

R         Thank you for clarifying some of the questions. Regarding the technical issues, I still have questions.



R         1)  I am obviously not part of the NMR community, but I still think there has to be an equation for the asymmetry parameter which can be understood by everybody. Do I understand correctly that you average the quantity given in (5) over all jumps from 2->3 during the simulation, but you give different weights depending on the initial (or final, or both?) position of the particle (whether it is close to an active site or not?)

A         There are no different weights assigned to particles in different pools. The quantity Mj is a concentration. In the simulations it is the number of particles in pool j. The quantity k_ij denotes the rate of transitions from pool j to pool i. The program counts the number of particles passing from pool  j to i and assigns that to k_ij M_j. The denominator of (5) is the total number of jumps from one pool to another including jumps within one pool (k_ii Mi). So, the asymmetry parameter defined in (5) is the number difference between forward and backward jumps divided by the total number of jumps. The total number of jumps is the sum of all differential jumps, i. e. the sum over all k_ij M_j. It is calulated in the program after completion of each simulation run and has been verified by comparison with the initially specified number of jumps.

R         2) Thanks for clarifying the details of the off-lattice simulation. It is of course ok and common practice to use simulation units and not SI units. Your simulation units are apparently defined in terms of the mass m of the particles, the particle diameter σ (I assume it is one in your units), and the energy (I assume you set kT=1 when setting up the Maxwell-Boltzmann velocity distribution). This defines the time unit τ=(m s² /kT)^1/2  . In these units, your time step is Δt = 1 τ, which is very large. In molecular dynamics simulations, typical values are Δt = 10^-3 τ or less, and having a too large time step can have a severe impact on the results.

A         We agree that the time step has been too large and as a result has introduced a false image of asymmetry by an unintentional overlap of the hard circles and by possibly skipping collisions. The simulation has now been tested at shorter time steps to determine if such an error was introduced. We find that the asymmetry parameter gets smaller as the time step is decreased. Furthermore, we are looking into the reviewer’s suggestion of modifying the gas diffusion code to accomplish event driven calculations.

R         On the other hand, you do not really describe a Molecular Dynamics simulation, it is rather another type of (off-lattice) Monte Carlo simulation. For example, your collisions preserve the energy of the two colliding particles, but not their momentum. As a Monte-Carlo simulation, it does not preserve detailed balance, and therefore, again, it is not surprising that the results also break detailed balance.

A         The gas phase simulation is a common time-driven elastic hard circle model with walls rather than periodic boundaries. An initial distribution of speeds is generated and the particles are given a random initial direction of travel. Immediately after a collision, new velocities and deflection angles for the two particles are determined from conservation of momentum and kinetic energy by:

As mentioned in the previous response, the asymmetry parameter decreases with shorter time step suggesting that detailed balance is recovering for the shorter times.  Including event driven dynamics into this code, as also mentioned above, should help in the study of detailed balance violation.

R         3) Regarding the comment "We found nonzero asymmetry parameters also when choosing the destination cell for a jump from all vacant neighbor cells at random." I would like to note that this algorithm also breaks detailed balance. In order to maintain detailed balance, you have to choose a jump randomly from all neighbor cells and then reject the move if the neighbor cell is filled. Rejecting means the particle stays where it is and does not move at all, but the follow-up configuration still counts for the overall statistics.  But I suppose the authors are aware of this, since they also state " Choosing an algorithm that satisfies detailed balance by construction precludes testing detailed balance" (a statement to which I fully agree.)

A         Thank you for clarifying! We reran the calculations for Fig. 6e choosing randomly from all neighbor positions, free and occupied. Indeed, the asymmetry parameter produces noise more than one order of magnitude lower than the values observed with zero probability assigned to randomly chosen jumps to occupied positions. We now understand, that by introducing a probability to jumps, detailed balance is violated. So, we are driving the imbalance by our vacancy diffusion algorithm.

R         4) Regarding the comment on currents: Cyclic exchange, as long as it is persistent and does not average to zero on the long run, would also count as current in the  definition of thermal equilibrium.

A         ok.

R         5) Regarding thermalization: This not only means thermalization with respect to velocities. Asking whether a system thermalizes is the same as asking whether a system reaches thermal equilibrium in the above sense, i.e., a stationary state without stationary (macroscopic or microscopic) currents.  Testing this is generally difficult in a simulation. The velocity distribution usually approaches the Maxwell-distribution very quickly, but other quantities usually equilibrate much more slowly. In your simulations, you could for example test the system without force terms and check whether the system ever reaches the Boltzmann distribution with respect to positions. In the absence of any forces, the particles should be uniformly distributed in the pore. This is probably not the case.

A         What happens in thermal equilibrium is at the heart of the question. Detailed balance demands no currents. Experimental evidence (Song 2000) and tested theory (Brownstein 1977) indicate the existence of diffusion eigenmodes of fluids confined to pores, which describe spatially oscillating distributions of nuclear magnetization components in thermal equilibrium and which decay in distinct ways under the impact of diffusion and the boundary geometry. Considering, that three-site exchange probes Fick’s second law, we interpret our observed currents, to be a pore-resonance effect on translational motion which obeys the diffusion equation.

Conclusion

Thanks to the lessons learned from the expert advice of the reviewer of our manuscript we understand that our manuscript needs major revision. The main point is that the reported computer simulations do not indicate a violation of the principle of detailed balance. Therefore, we did not find any indication that three-site NMR exchange maps are asymmetric in thermodynamic equilibrium, and the observed asymmetry needs to be attributed to experimental or data processing (e.g. inverse Laplace transformation) deficiencies. A second outcome of the computer simulations remains without major revision. This is that the observed particle dynamics obey the diffusion equation and are linked to diffusion eigenmodes with the consequence that diffusion eigenmodes may possibly be driven by external stimuli like the violin bow enforcing resonance vibration of a Chladni plate. This is a technically interesting perspective as heterogeneous catalysis may possibly be enhanced be oscillating electrical, magnetic, or mechanical (ultrasound) fields. Pending the editor’s decision, we would revise our manuscript accordingly.

Thank you for your persistence and for taking the time to educate us on computer modelling. Obviously, our expertise is not in this area. Nevertheless, we hoped to gain some insight into the symmetry properties of relaxation-exchange maps in NMR in this way.

R         Thank you for clarifying some of the questions. Regarding the technical issues, I still have questions.



R         1)  I am obviously not part of the NMR community, but I still think there has to be an equation for the asymmetry parameter which can be understood by everybody. Do I understand correctly that you average the quantity given in (5) over all jumps from 2->3 during the simulation, but you give different weights depending on the initial (or final, or both?) position of the particle (whether it is close to an active site or not?)

A         There are no different weights assigned to particles in different pools. The quantity Mj is a concentration. In the simulations it is the number of particles in pool j. The quantity k_ij denotes the rate of transitions from pool j to pool i. The program counts the number of particles passing from pool  j to i and assigns that to k_ij M_j. The denominator of (5) is the total number of jumps from one pool to another including jumps within one pool (k_ii Mi). So, the asymmetry parameter defined in (5) is the number difference between forward and backward jumps divided by the total number of jumps. The total number of jumps is the sum of all differential jumps, i. e. the sum over all k_ij M_j. It is calulated in the program after completion of each simulation run and has been verified by comparison with the initially specified number of jumps.

R         2) Thanks for clarifying the details of the off-lattice simulation. It is of course ok and common practice to use simulation units and not SI units. Your simulation units are apparently defined in terms of the mass m of the particles, the particle diameter σ (I assume it is one in your units), and the energy (I assume you set kT=1 when setting up the Maxwell-Boltzmann velocity distribution). This defines the time unit τ=(m s² /kT)^1/2  . In these units, your time step is Δt = 1 τ, which is very large. In molecular dynamics simulations, typical values are Δt = 10^-3 τ or less, and having a too large time step can have a severe impact on the results.

A         We agree that the time step has been too large and as a result has introduced a false image of asymmetry by an unintentional overlap of the hard circles and by possibly skipping collisions. The simulation has now been tested at shorter time steps to determine if such an error was introduced. We find that the asymmetry parameter gets smaller as the time step is decreased. Furthermore, we are looking into the reviewer’s suggestion of modifying the gas diffusion code to accomplish event driven calculations.

R         On the other hand, you do not really describe a Molecular Dynamics simulation, it is rather another type of (off-lattice) Monte Carlo simulation. For example, your collisions preserve the energy of the two colliding particles, but not their momentum. As a Monte-Carlo simulation, it does not preserve detailed balance, and therefore, again, it is not surprising that the results also break detailed balance.

A         The gas phase simulation is a common time-driven elastic hard circle model with walls rather than periodic boundaries. An initial distribution of speeds is generated and the particles are given a random initial direction of travel. Immediately after a collision, new velocities and deflection angles for the two particles are determined from conservation of momentum and kinetic energy by:

As mentioned in the previous response, the asymmetry parameter decreases with shorter time step suggesting that detailed balance is recovering for the shorter times.  Including event driven dynamics into this code, as also mentioned above, should help in the study of detailed balance violation.

R         3) Regarding the comment "We found nonzero asymmetry parameters also when choosing the destination cell for a jump from all vacant neighbor cells at random." I would like to note that this algorithm also breaks detailed balance. In order to maintain detailed balance, you have to choose a jump randomly from all neighbor cells and then reject the move if the neighbor cell is filled. Rejecting means the particle stays where it is and does not move at all, but the follow-up configuration still counts for the overall statistics.  But I suppose the authors are aware of this, since they also state " Choosing an algorithm that satisfies detailed balance by construction precludes testing detailed balance" (a statement to which I fully agree.)

A         Thank you for clarifying! We reran the calculations for Fig. 6e choosing randomly from all neighbor positions, free and occupied. Indeed, the asymmetry parameter produces noise more than one order of magnitude lower than the values observed with zero probability assigned to randomly chosen jumps to occupied positions. We now understand, that by introducing a probability to jumps, detailed balance is violated. So, we are driving the imbalance by our vacancy diffusion algorithm.

R         4) Regarding the comment on currents: Cyclic exchange, as long as it is persistent and does not average to zero on the long run, would also count as current in the  definition of thermal equilibrium.

A         ok.

R         5) Regarding thermalization: This not only means thermalization with respect to velocities. Asking whether a system thermalizes is the same as asking whether a system reaches thermal equilibrium in the above sense, i.e., a stationary state without stationary (macroscopic or microscopic) currents.  Testing this is generally difficult in a simulation. The velocity distribution usually approaches the Maxwell-distribution very quickly, but other quantities usually equilibrate much more slowly. In your simulations, you could for example test the system without force terms and check whether the system ever reaches the Boltzmann distribution with respect to positions. In the absence of any forces, the particles should be uniformly distributed in the pore. This is probably not the case.

A         What happens in thermal equilibrium is at the heart of the question. Detailed balance demands no currents. Experimental evidence (Song 2000) and tested theory (Brownstein 1977) indicate the existence of diffusion eigenmodes of fluids confined to pores, which describe spatially oscillating distributions of nuclear magnetization components in thermal equilibrium and which decay in distinct ways under the impact of diffusion and the boundary geometry. Considering, that three-site exchange probes Fick’s second law, we interpret our observed currents, to be a pore-resonance effect on translational motion which obeys the diffusion equation.

Conclusion

Thanks to the lessons learned from the expert advice of the reviewer of our manuscript we understand that our manuscript needs major revision. The main point is that the reported computer simulations do not indicate a violation of the principle of detailed balance. Therefore, we did not find any indication that three-site NMR exchange maps are asymmetric in thermodynamic equilibrium, and the observed asymmetry needs to be attributed to experimental or data processing (e.g. inverse Laplace transformation) deficiencies. A second outcome of the computer simulations remains without major revision. This is that the observed particle dynamics obey the diffusion equation and are linked to diffusion eigenmodes with the consequence that diffusion eigenmodes may possibly be driven by external stimuli like the violin bow enforcing resonance vibration of a Chladni plate. This is a technically interesting perspective as heterogeneous catalysis may possibly be enhanced be oscillating electrical, magnetic, or mechanical (ultrasound) fields. Pending the editor’s decision, we would revise our manuscript accordingly.

Thank you for your persistence and for taking the time to educate us on computer modelling. Obviously, our expertise is not in this area. Nevertheless, we hoped to gain some insight into the symmetry properties of relaxation-exchange maps in NMR in this way.

R         Thank you for clarifying some of the questions. Regarding the technical issues, I still have questions.



R         1)  I am obviously not part of the NMR community, but I still think there has to be an equation for the asymmetry parameter which can be understood by everybody. Do I understand correctly that you average the quantity given in (5) over all jumps from 2->3 during the simulation, but you give different weights depending on the initial (or final, or both?) position of the particle (whether it is close to an active site or not?)

A         There are no different weights assigned to particles in different pools. The quantity Mj is a concentration. In the simulations it is the number of particles in pool j. The quantity k_ij denotes the rate of transitions from pool j to pool i. The program counts the number of particles passing from pool  j to i and assigns that to k_ij M_j. The denominator of (5) is the total number of jumps from one pool to another including jumps within one pool (k_ii Mi). So, the asymmetry parameter defined in (5) is the number difference between forward and backward jumps divided by the total number of jumps. The total number of jumps is the sum of all differential jumps, i. e. the sum over all k_ij M_j. It is calulated in the program after completion of each simulation run and has been verified by comparison with the initially specified number of jumps.

R         2) Thanks for clarifying the details of the off-lattice simulation. It is of course ok and common practice to use simulation units and not SI units. Your simulation units are apparently defined in terms of the mass m of the particles, the particle diameter σ (I assume it is one in your units), and the energy (I assume you set kT=1 when setting up the Maxwell-Boltzmann velocity distribution). This defines the time unit τ=(m s² /kT)^1/2  . In these units, your time step is Δt = 1 τ, which is very large. In molecular dynamics simulations, typical values are Δt = 10^-3 τ or less, and having a too large time step can have a severe impact on the results.

A         We agree that the time step has been too large and as a result has introduced a false image of asymmetry by an unintentional overlap of the hard circles and by possibly skipping collisions. The simulation has now been tested at shorter time steps to determine if such an error was introduced. We find that the asymmetry parameter gets smaller as the time step is decreased. Furthermore, we are looking into the reviewer’s suggestion of modifying the gas diffusion code to accomplish event driven calculations.

R         On the other hand, you do not really describe a Molecular Dynamics simulation, it is rather another type of (off-lattice) Monte Carlo simulation. For example, your collisions preserve the energy of the two colliding particles, but not their momentum. As a Monte-Carlo simulation, it does not preserve detailed balance, and therefore, again, it is not surprising that the results also break detailed balance.

A         The gas phase simulation is a common time-driven elastic hard circle model with walls rather than periodic boundaries. An initial distribution of speeds is generated and the particles are given a random initial direction of travel. Immediately after a collision, new velocities and deflection angles for the two particles are determined from conservation of momentum and kinetic energy by:

As mentioned in the previous response, the asymmetry parameter decreases with shorter time step suggesting that detailed balance is recovering for the shorter times.  Including event driven dynamics into this code, as also mentioned above, should help in the study of detailed balance violation.

R         3) Regarding the comment "We found nonzero asymmetry parameters also when choosing the destination cell for a jump from all vacant neighbor cells at random." I would like to note that this algorithm also breaks detailed balance. In order to maintain detailed balance, you have to choose a jump randomly from all neighbor cells and then reject the move if the neighbor cell is filled. Rejecting means the particle stays where it is and does not move at all, but the follow-up configuration still counts for the overall statistics.  But I suppose the authors are aware of this, since they also state " Choosing an algorithm that satisfies detailed balance by construction precludes testing detailed balance" (a statement to which I fully agree.)

A         Thank you for clarifying! We reran the calculations for Fig. 6e choosing randomly from all neighbor positions, free and occupied. Indeed, the asymmetry parameter produces noise more than one order of magnitude lower than the values observed with zero probability assigned to randomly chosen jumps to occupied positions. We now understand, that by introducing a probability to jumps, detailed balance is violated. So, we are driving the imbalance by our vacancy diffusion algorithm.

R         4) Regarding the comment on currents: Cyclic exchange, as long as it is persistent and does not average to zero on the long run, would also count as current in the  definition of thermal equilibrium.

A         ok.

R         5) Regarding thermalization: This not only means thermalization with respect to velocities. Asking whether a system thermalizes is the same as asking whether a system reaches thermal equilibrium in the above sense, i.e., a stationary state without stationary (macroscopic or microscopic) currents.  Testing this is generally difficult in a simulation. The velocity distribution usually approaches the Maxwell-distribution very quickly, but other quantities usually equilibrate much more slowly. In your simulations, you could for example test the system without force terms and check whether the system ever reaches the Boltzmann distribution with respect to positions. In the absence of any forces, the particles should be uniformly distributed in the pore. This is probably not the case.

A         What happens in thermal equilibrium is at the heart of the question. Detailed balance demands no currents. Experimental evidence (Song 2000) and tested theory (Brownstein 1977) indicate the existence of diffusion eigenmodes of fluids confined to pores, which describe spatially oscillating distributions of nuclear magnetization components in thermal equilibrium and which decay in distinct ways under the impact of diffusion and the boundary geometry. Considering, that three-site exchange probes Fick’s second law, we interpret our observed currents, to be a pore-resonance effect on translational motion which obeys the diffusion equation.

Conclusion

Thanks to the lessons learned from the expert advice of the reviewer of our manuscript we understand that our manuscript needs major revision. The main point is that the reported computer simulations do not indicate a violation of the principle of detailed balance. Therefore, we did not find any indication that three-site NMR exchange maps are asymmetric in thermodynamic equilibrium, and the observed asymmetry needs to be attributed to experimental or data processing (e.g. inverse Laplace transformation) deficiencies. A second outcome of the computer simulations remains without major revision. This is that the observed particle dynamics obey the diffusion equation and are linked to diffusion eigenmodes with the consequence that diffusion eigenmodes may possibly be driven by external stimuli like the violin bow enforcing resonance vibration of a Chladni plate. This is a technically interesting perspective as heterogeneous catalysis may possibly be enhanced be oscillating electrical, magnetic, or mechanical (ultrasound) fields. Pending the editor’s decision, we would revise our manuscript accordingly.