Given two Hamiltonians \(H_0,H_1=H_0+V\) corresponding respectively to a free particle Hamiltonian and the free particle Hamiltonian plus its interaction energy with a scattering centre, we can compute the S-matrix, i.e., the matrix that determines the quantum mechanical amplitude of the process involving a particle coming from a free particle state at an infinite distance from the scattering centre at time \(t=-\infty\), entering into an initial scattered state, interacting with the scattering centre, entering into a final scattered state, and then escaping away to infinity at time \(t\rightarrow\infty\) to a final free particle state. Such amplitudes were first studied by Lippmann and Schwinger (Steven Weinberg, The Quantum Theory of Fields, Vol. 1, W.O. Amrein, Hilbert Space Methods in Quantum Mechanics) and then analyzed in a mathematically precise way by [T. Kato], [W. Amrein], and [K.B. Sinha]. In this paper, we first derive the Lippmann-Schwinger equation for the input and output scattered states in terms of the input and output free particle states using intuitive arguments originally due to Lippmann and Schwinger, and then, using rigorous operator theoretic arguments revolving around the spectral theorem for unbounded self-adjoint operators in a Hilbert space, we obtain an elegant formula for the scattering matrix elements in terms of the interaction potential and the free particle Hamiltonian between two free particle states corresponding to the same measure. We then discuss a computationally more efficient method based on the Dyson series expansion (Steven Weinberg, The Quantum Theory of Fields, Vol. 1), which works even in the case when the interaction Hamiltonian is time varying, in contrast to the Lippmann-Schwinger method, for calculating the S-matrix. We assume that the scattering potential, which can be time varying, can be controlled by incorporating control parameters and explain how the resulting unitary S-matrix can be made to approximate a given unitary gate in infinite-dimensional Hilbert space by optimizing the Hilbert-Schmidt/Frobenius norm distance between the given gate and the S-matrix gate, in which the latter is approximated by a truncated Dyson series. Of course, the S-matrix given by the Dyson series is a Taylor series functional expansion in the interaction potential and hence is a highly nonlinear function of the control parameters, and hence we propose numerical methods for such optimization. We then observe that if the control parameters are taken as random variables with joint moments, then the TPCP map involving statistical averaging of the adjoint action of the S-matrix acting on an initial mixed state w.r.t. the probability distribution of the control parameters will result in a final mixed state that is linear in the statistical moments of the control parameters, and hence optimization w.r.t. these moments can be easily carried out. Our method of TPCP map design is to optimize the sum of Frobenius distance squares between a sequence of desired output mixed states and the S-matrix based formula for the TPCP map acting on the corresponding input mixed states w.r.t. the statistical moments of the control parameters. We then generalize this method to include the case of TPCP maps obtained from the Hudson-Parthasarathy quantum noisy Schrödinger equation. In summary, we simulate the S-matrix gate for independent realizations of the control parameters having a probability distribution defined by their joint moments obtained by the above optimization procedure and apply the adjoint action of these independent S-matrices on a given input state, and take their ensemble average to obtain a good approximation to the action of the desired TPCP map on any given input state. This suggests that TPCP maps of arbitrarily large size can be realized using scattering experiments in the laboratory.