Aeroelastic response at bumblebee scale is explored using a well-validated Navier-Stokes equation solver, fully-coupled with a structural dynamics solver. Hover flight at Re=1.0×10 3 is considered, which results in a purely passive wing rotation due to the dynamic balance between aerodynamic loading, elastic restoring force, and inertial force of the wing. A systematic study with effects of variation of motion amplitude and wing stiffness on the resulting aerodynamic and structural dynamic response is reported. The aeroelastic response is non-periodic and varies cycle to cycle. Highest time-averaged lift of 1.43 is obtained at a moderate amplitude and lowest wing stiffness. Optimal efficiency corresponds to the largest motion amplitude with a moderate wing stiffness. Dual vortical structures are observed at both ends of the wing during a stroke, leading to multiple wing-wake interactions. The vortical evolution is more chaotic than at Re=1.0×10 2 that corresponds to the fruit fly scale. Nevertheless, the time history of lift for optimal efficiency motions are remarkably similar for both scales. Moreover, the time averaged lift scales with the shape deformation parameter, reinforcing the observation made at both fruit fly and water tunnel scales. Finally, the reduced frequency for the optimal efficiency motion is close to experimentally observed values for hovering bumblebees, suggesting that bumblebee wing kinematics may aim to be aerodynamically optimal. Nomenclature c = Chord [m] C L = Coefficient of lift [1] C P = Coefficient of power input [1] E = Young's Modulus [Pa] f = Motion frequency [1/s] f 1 = First natural frequency of the wing [1/s] F = Fluid force acting on the wing per unit length [N/m] h = Plunge motion of the wing [m] h a = Plunge amplitude [m] h s = Thickness of the wing [m] k = Reduced frequency, / fc U [1] p = Pressure [Pa] Re = Reynolds number, / Uc [1] St = Strouhal number, ~/ a fh U [1] t = Time [s] U = Maximum plunge velocity: 2 a fh for hover [m/s 2 ] u = Velocity [m/s] v = Wing displacement: v=w+h [m] w = Wing displacement relative to the imposed motion h [m] = Passive pitch angle [degrees] a = Passive pitch amplitude [degrees] m = Mid-stroke passive pitch angle [degrees] e = End-of-the-stroke passive pitch angle [degrees] Φ = Stroke angle [degrees] = Phase lag between passive pitch and plunge motion [degrees] = Non-dimensional relative shape parameter [1] = Propulsive efficiency: / L P C C [1] = Kinematic viscosity of the fluid [m 2 /s] 0 = Effective inertia: * * 2 ( / ) s h k [1] 1 = Effective stiffness: * 2 / (12 ) s f Eh U [1] f = Density of the fluid [kg/m 3 ] s = Density of the wing structure [kg/m 3 ] = Vorticity, non-zero component of u [1/s] () * = Non-dimensional variables ( ) = Time-averaged variables